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ABSTRACT 


A higher-order composite box beam theory is developed to model rectangular beams 
with arbitrary wall thicknesses. The theory, which is based on a refined displacement 
field, is a three-dimensional model which approximates the elasticity solution so that the 
beam cross-sectional properties are not reduced to one-dimensional beam parameters. Both 
inplane and out-of-plane warping are included automatically in the formulation. The model 
can accurately capture the transverse shear stresses through the thickness of each wall while 
satisfying stress free boundary conditions on the inner and outer surfaces of the beam. 

Numerical results are presented for beams with different wall thicknesses and aspect 
ratios. The static results are correlated with available experimental data and show excellent 
agreement. The dynamic results which are correlated with a general purpose finite element 
code show the importance of including inplane and out-of-plane warping deformations in 
the formulation. 

The developed beam theory is then used to model the load carrying member of a tilt- 
rotor blade which has thick-walled sections. A procedure is developed for computing the 
aeroelastic stability of the tilt-rotor blade based on the composite box beam model. The 
aerodynamic loads are calculated using the classical blade element momentum theory. 
Analytical expressions for the lift and drag are obtained based on the blade planform with 
corrections for the high lift capability of rotor blades. The aerodynamic analysis is coupled 
with the structural model to formulate the complete coupled equations of motion for 
aeroelastic analyses. 

Finally, a multidisciplinary optimization procedure is developed to improve the 
aerodynamic, structural and aeroelastic performance of the tilt-rotor aircraft. The objective 
functions include the figure of merit in hover and the high speed cruise propulsive 
efficiency. Structural, aerodynamic and aeroelastic stability criteria are imposed as 
constraints on the problem. The Kreisselmeier-Steinhauser function is used to formulate 

iii 



the multiobjective function problem and a hybrid approximate analysis is used to reduced 
the computational effort. The search direction is determined by the Broyden-Fletcher- 
Goldfarb-Shanno algorithm. The optimum results are compared with the baseline values 
and show significant improvements in the overall performance of the tilt-rotor blade. 
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1. Introduction 


The tilt-rotor aircraft has become an increasingly viable concept over the last several 
years [1-4]. The goal is to combine the hover performance of a helicopter in take-off and 
landing with fixed-wing aircraft like efficiencies in high speed cruise. The aircraft is 
similar in appearance to conventional fixed-wing aircraft, however, large diameter rotors 
are tip mounted on the wings. These rotors are mounted on pylon assemblies which are 
capable of rotation through 90 degrees so that the aircraft can convert between the various 
flight regimes (e.g. hover, transition and cruise) required of this aircraft (Figs. 1 - 3). 



Fig. 1.1 XV- 15 tilt-rotor in helicopter mode. 



Fig. 1.2 XV- 15 tilt-rotor in transition/conversion mode. 


There are many design requirements associated with tilt-rotor performance including a 
low disc loading in the hover configuration and the ability to rotate the rotors forward to 
achieve cruise speeds up to 450 knots [5]. The problem becomes more complex since in 
vertical flight and in hover, a large portion of the rotor is directly over the wing which 


2 


produces a large downwash effect upon the wing. The down wash effect increases thrust 
requirement of the aircraft by approximately 10-12 percent [6]. Other problems associated 
with this configuration are related to high helical tip Mach numbers (M t i p ) which represent 
a critical performance issue in high speed cruise (350 - 450 knots). Aeroelastic stability is 
an important consideration in the design of tilt-rotors. Due to the large thrust requirement 
in hover, the prop-rotors have a much greater radius than standard propellers. This 
increases the tip speed which in cruise may cause individual blade flutter or a coupled 
flexible motion between the rotor, wing and pylon known as whirl flutter. Since civil tilt- 
rotors are required to be stable at a 20 percent margin above their dive speed (defined to be 
15 percent larger than the cruise speed), this means that the flutter speed must be larger than 
621 knots for a tilt-rotor with a cmise speed of 450 knots. 



There are several different techniques which can be used to address these issues. For 
example, the tip Mach number can be reduced through rotor tip speed reduction or through 
the use of blade sweep which reduces the effective Mach number. Another alternative is to 
increase the drag divergence Mach number (Mdd) at the tip to values above Mtip. This can 
be accomplished through reductions in the blade thickness. However, each of these 
options will adversely affect the hover performance, drive system weight or aeroelastic 
stability of the rotor blade. In the helicopter mode, to maintain a high figure of merit in 
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hover, the solidity of the blade must be increased since thinner airfoils are used for 
maintaining efficiency in cruise. 

1.1 Tilt-Rotor Design 

Due to the many conflicting requirements imposed on prop-rotor performance between 
hover, conversion and airplane mode, the use of formal numerical optimization techniques 
is well suited for studying the design trade-offs associated with such aircraft. Although 
there has been a significant amount of research performed on the optimization of rotary 
wing aircraft, only very few studies have investigated tilt-rotor design issues. Recently, 
research efforts have been initiated by Chattopadhyay et al. [7-16] to develop formal 
optimization techniques to address these issues. In Refs. 7-9, optimization procedures 
were developed to maximize the high speed cruise propulsive efficiency without degrading 
the hover figure of merit. An optimization procedure was developed in Ref. 10 to address 
the problem of aeroelastic stability of high speed proprotors. In Ref. 1 1, the drive system 
weight was minimized and the associated trade-offs in cruise efficiency were investigated. 
The integrated aerodynamic, aeroelastic and structural optimization problem was addressed 
in Ref. 12. A purely aerodynamic multiobjective optimization procedure for improved high 
speed cruise and hovering performance using planform and airfoil characteristics as design 
variables was developed by McCarthy et al. [13]. In Refs. 14 and 15, a two level 
decomposition technique was developed by Chattopadhyay et al. for the combined 
aerodvnamic/structural design of high speed prop-rotor blades. At the upper level, the 
aerodynamic performance was improved using continuous optimization techniques. The 
structural performance was improved at the lower level using composite ply orientations a 
of thin-walled box beam as design variables. More recently, McCarthy and Chattopadhyay 
[16], developed a procedure for a combined wing/rotor optimization study of high speed 
tilt-rotors. The analysis was based on comprehensive rotorcraft algorithms and design 
criteria from both the rotor and the wing were included in the formulation. 
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1.2 Aeroelastic Stability 

Rotor dynamics and aeroelastic response are critical to the design of vertical take-off 
and landing and short take-off (V/STOL) aircraft. In 1966, Reed [17] described the basic 
phenomenon of propeller whirl flutter instability for V/STOL aircraft. Since this time, 
several survey papers have attempted to address the critical issues associated with 
aeroelastic stability. In 1976, Kvaternik [18] noted that the blade inplane flexibility can 
have a significant effect on the stability of the system. More recently, Friedmann [19] 
described the necessary requirements in aeroelastic stability modeling. The significant 
potential for aeroelastic tailoring from using composite rotors in a comprehensive 
aeroelastic stability analysis was discussed. 

There has also been a significant amount of in depth research efforts to investigate the 
physics of aeroelastic stability. Edenborough [20] investigated the stability of rotor-pylon 
configuration for a tilt-rotor aircraft. It was concluded from the study that stable 
rotor/pylon configuration can be designed for high speed operation. The configuration 
studied in the paper was modeled analytically up to speeds of 250 knots which was then 
validated by experimental data up to a speed of about 195 knots. Kaza [21] investigated the 
effect of steady-state coning angle and damping of the flapping hinge of the blades on the 
whirl flutter stability boundary. This analytical study was based on one-dimensional beam 
analysis. More recently, Johnson [22,23] has identified the modeling requirements to 
accurately predict aeroelastic stability. From these studies it was concluded that an accurate 
structural representation of the blade is essential to properly model dynamic stability. 

Johnson [24] calculated the performance, loads and stability of the XV- 15 tilt-rotor and 
compared the results to wind-tunnel and flight test measurements to assess the requirements 
for additional experimental data and further analytical development. The basic dynamic 
problems of advanced rotor system, including aeroelasticity of tilt-rotor blades, was also 
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investigated by Johnson in Ref. 25. In this study the need for accurate elastic 
representation of the blade dynamics was established. 

1.3 Aeroelastic Tailoring 

Due to the fact that aeroelastic stability is critical in the design of rotary wing aircraft, a 
significant amount of research has been performed investigating aeroelastic tailoring for 
rotary wing aircraft over the years. The aeroelastic tailoring of composite rotor blades 
based on the Classical Laminate Theory (CLT) for a thin-walled cylindrical tube was 
investigated by Mansfield and Sobey [26]. A detailed investigation of aeroelastic stability 
of helicopter blades was performed by Chopra and Sivaneri [27]. Quasi- steady airloads 
were used and the equations of motion were linearized using a small perturbation theory. 
The equations of motion were then solved using a finite element technique based on one- 
dimensional isotropic beam approach. 

The aeroelastic stability of composite helicopter blades was investigated using finite 
element techniques by Hong and Chopra [28]. The structural element in the blade was 
modeled as a thin-walled composite box beam. Only axial and inplane shear stresses were 
included in the formulation. This approach was then modified by the authors to investigate 
the stability of a bearingless rotor by modeling the rotor flex beam as an I-beam [29]. 

Rosen and Rand [30] developed a theory which describes the aeroelastic behavior of 
curved helicopter blades for hovering and axial flight. In this approach, the blade was 
modeled as a solid, isotropic beam. The control of tilt-rotor flutter was analyzed by Nasu 
[31] using a rotor model consisting of a straight, fixed wing, a pylon attached to the wing 
tip and a three-bladed rotor. Each rotor blade had two bending degrees of freedom and one 
torsional degree of freedom about the elastic axis. Torok and Chopra [32] used a 
comprehensive rotor aeroelastic analysis to investigate the stability of a hingeless helicopter 
blade. The aeroelastic response and blade loads of a helicopter blade were investigated by 
Smith and Chopra [33]. A thin-walled box beam model based on the CLT was used as the 
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load carrying member of the blade in this study. Using the same structural model, Smith 
[34] performed a parametric study to investigate composite tailoring effects on helicopter 
blade vibration and flutter. 

The influence of several system design parameters on tilt-rotor aeroelastic stability in 
high speed cruise was addressed by Nixon using a parametric study [35]. The study 
indicates that the separation of the beam and torsion frequencies can have a significant 
effect on the stability of the beam mode. It was also shown that the rotor thrust level has a 
negligible effect on the stability. This study was later enhanced by the author to investigate 
possibilities of elastic tailoring using composite rotors [36]. 

A investigation of aeroelastic stability of an advanced geometry helicopter blade which 
includes fuselage dynamic effects included in the formulation was performed by Bir and 
Chopra [37]. The aeromechanical stability of bearingless composite helicopter blade in 
forward flight was investigated by Tracy and Chopra [38]. Transverse shear effects were 
included in the formulation although it was assumed that the transverse displacements can 
be decoupled into bending terms and shear terms. A further assumption made in the study 
is that the transverse shear forces were equal to zero. 

Nonlinear, large amplitude aeroelastic behavior of composite rotor blades was 
investigated using one-dimensional beam theory by Kim and Dugundji [39]. The nonlinear 
equations of motion were solved using a Newton-Raphson technique. Yen [40] 
investigated the effects of the blade tip shape on rotor dynamics and aeroelastic response 
using a parametric study. Several conflicting requirements were observed including the 
fact that the blade weight had to be increased in order to achieve reductions in the vibratory 
behavior. The aeroelastic analysis of a fixed wing and a rotating wing in hover were 
addressed by Nibbeling and Peters [41] using inflow dynamics and a linear structural 
model which included only elastic bending and torsion. The main focus of this study was 
inflow dynamics and the structural model was relatively simplistic in nature. Yuan and 
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Friedmann [42] used an anisotropic one-dimensional beam model to investigate aeroelastic 
stability of swept tip composite helicopter blades. A detailed finite element model was 
required to reduce the three-dimensional behavior of the blade to one-dimensional beam like 
parameters. Only extensional and transverse shear stresses were included in the 
formulation. 

The use of composite materials also allows tailoring capabilities which can be used to 
improve the rotor structural and dynamic behavior. Structural tailoring using composite 
blades has recently been investigated as a means for improving the aeroelastic stability of 
rotary wing aircraft [33,35,36,42-44]. Parametric studies have been performed to 
investigate composite tailoring in rotary wing applications in Refs. 33, 35, 36 and 42. Of 
these, only Refs. 35 and 36 investigate tilt-rotor aircraft. Formal optimization techniques 
have been developed by Ganguli and Chopra for the optimization of composite rotor blades 
in helicopters [43,44]. In these studies, continuous values for the ply angles were used as 
design variables. 

Improvements in aeroelastic stability through structural tailoring without consideration 
of other requirements such as aerodynamic and structural performance can lead to 
unrealistic designs. These criteria can be effectively addressed using multidisciplinary 
optimization techniques. 

1.4 Structural Modeling 

For the tilt-rotor, the combination of the low disk-loading of the prop-rotors and the 
high inflow ratio (the ratio of axial velocity to rotor tip speed) make the dynamics and 
aeroelastic response of this configuration unique and more complex than that of a 
helicopter. It is therefore very important that the issues of individual blade flutter and whirl 
flutter be addressed very early in the design process. As noted in the literature, this can be 
accomplished by structurally tailoring the rotor for proper stiffness and elastic couplings in 
an effort to raise blade frequencies and minimize inplane motions at the rotor hub. 
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However, in order to adequately model the dynamics and aeroelasticity of tilt-rotors, it is 
essential that an accurate structural model of the rotor be incorporated into the analysis 
[23,24]. Since efficient structural modeling of the rotor is key to the overall analysis of the 
aircraft, research efforts have been initiated to investigate this issue. 

1.4.1 Beam theories : Beam theories associated with isotropic materials have been well 
understood for years and these theories tend to predict the structural and dynamic response 
quite accurately [45,46]. In recent years, the use of composite materials in rotary wing 
applications has increased due to the favorable strength-to-weight characteristics of these 
anisotropic materials. However, beam theories for anisotropic materials are not as well 
established as they are for isotropic materials. Recently some research has been reported in 
deriving composite beam theories [47-80]. 

A complex finite element methodology capable of modeling anisotropic beams with 
arbitrary cross sections was developed by Womdle [47]. The method was based on a fully 
three-dimensional finite element theory in which the cross-sectional warping was 
superimposed on the assumed displacement field. A comprehensive theory for modeling 
anisotropic material was developed by Giavotto et al. [48]. In this methodology, the three- 
dimensional behavior of the beam was reduced to one-dimensional beam like parameters 
based on a cross-sectional analysis. St. Venant type warping was superimposed on the 
assumed displacement field. 

Krenk and Gunneskov [49] formulated a theory for pre-twisted turbine blades which 
included finite shear flexibility. The pre-twist was accounted for, in this theory, through 
the axial derivative of the St. Venant warping function and the shear stresses were 
decoupled into torsion and shear contributions. Using this approach, the shear flexibility 
could only be approximated in blades with moderate wall thicknesses. 

Another beam theory for anisotropic materials was due to Bauchau [50]. Out-of-plane 
warping was included in the formulation but the model assumed that each section was 
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infinitely rigid in its own plane, therefore inplane warping was neglected. The out-of-plane 
warping was based on so-called “eigenwarping functions” which was an improvement to 
the St. Venant warping. Only the axial strain and torsional related shear strain were 
included in the formulation. The theory was extended by the author to include the large 
displacement analysis of naturally curved and twisted composite beams [51]. Although 
shearing deformations and torsion related warping were included in this formulation, it was 
assumed that the cross sections do not deform in their own plane. Further, the beam theory 
was only derived for thin- walled beams. 

The static and dynamic behavior of helicopter blades using a finite element approach 
was addressed by Bauchau and Hong [52]. This study represented a first step towards 
developing a complete aeroelastic analysis. However, the theory used in this report was 
based on the Classical Laminate Theory which is valid only for thin-walled beams. 
Bauchau and Hong later developed a nonlinear composite beam theory to model helicopter 
rotor blades [53]. The modeling included arbitrarily large displacements and rotations but 
assumed small strains. Although warping effects were included, only the axial strain and 
the inplane and transverse shear strain appeared in the formulation. Shearing and warping 
deformations were investigated by Bauchau et al. [54] for a thin-walled composite 
sandwich beam and were compared with experimental data. The study once again assumed 
thin- walled sections. 

A nonlinear analysis of pre-twisted rods was developed by Rosen et al. [55]. Small 
strains were assumed in the formulation. This theory was later extend for modeling the 
dynamics of moving and rotating rods [56-57]. Large deformations were included in the 
formulation. A finite element model for the analysis of composite box beams was 
developed by Stemple [58]. The warping function, in this formulation, was assumed to be 
only out-of-plane and it was superimposed upon a displacement field consisting of three 
translations and three rotations. The warping was determined based upon a two- 
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dimensional cross-sectional analysis. The results from this theory have been shown to be 
as accurate as a full three-dimensional solution for thin-walled box beams [59]. 

Standard means of representing finite rotation in rigid-body kinematics, including 
orientation angles, Euler parameter and Rodrigues parameters were reviewed and compared 
by Hodges [60]. General kinematical relations for a beam which include moderate 
rotations were presented. Hinnant and Hodges [61] used a finite element methodology, 
which included geometric nonlinearities, to obtain results for a cantilevered beam with a tip 
mass which were compared to experimental values. A nonlinear formulation for the 
dynamics of initially curved and twisted beams in a moving frame was presented by 
Hodges [62], Both inplane and out-of-plane St. Venant warping displacements were 
assumed. The three-dimensional beam behavior was reduced to one-dimensional beam 
parameters and only unrestrained warping effects were included in the formulation. 
Rehfield et al. [63] investigated the nonclassical behavior of thin-walled composite beams 
with closed cross sections. The nonclassical behavior refers to elastic bending-shear 
coupling and restrained torsional warping. A decay length parameter was defined which 
approximates the effects of the restrained warping. 

Several refined composite beam theories based on the Variational Asymptotical method 
[64] have been developed by Hodges et al. [65-70]. The variational asymptotical method is 
a mathematical technique by which the three-dimensional analysis of the composite beam 
deformation can be split into a linear, two-dimensional cross-sectional analysis and a 
nonlinear, one -dimensional analysis. An ordering scheme is required to identify “small” 
terms which are then eliminated from the strain energy formulation. The method requires 
the energy functional to be expanded in terms of a small parameter and the theory is 
therefore truly only valid for thin- walled beams. 

An analytical model was developed by Kosmatka [71] for assessing the extension- 
bending-torsion coupling effects associated with anisotropic beams having non 
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homogeneous irregular cross sections and initial twist. The model was based on beam 
theory which reduces the three-dimensional behavior to one-dimensional beam equations 
based on a cross-sectional analysis. The vibration analysis of composite turbopropellers 
using a nonlinear beam-type finite element approach was studied by Kosmatka and 
Friedmann [72]. Constitutive relations were obtained by setting the three stresses within 
the cross section equal to zero, thus assuming that the cross section is rigid in its own 
plane. In Ref. 73, an analytical model was presented by Kosmatka and Dong for 
determining the displacement and stress distributions of the St. Venant extension, bending, 
torsion and flexural problem for prismatic, anisotropic beams of arbitrary cross sections. 
The problem was reduced to a state of plane stress. However, using this approach a finite 
element analysis must be performed for each unique cross section before the beam 
equations of motion can be solved. The behavior of a tip-loaded cantilever beam with an 
arbitrary cross section was investigated by Kosmatka using a power series solution for the 
out-of-plane flexural and torsional warping [74,75]. For complex cross sections, the 
warping results represented a best fit approximation to the exact St. Venant warping 
function. A beam theory based on a first-order displacement field with superimposed 
warping functions was used in the formulation which requires a complex two-dimensional 
finite element analysis of the cross section. 

A theoretical modeling of slender composite rotating beams was developed by Rand 
[76,77]. In addition to the classical degrees of freedom used in this theory, a three- 
dimensional warping field was superimposed on the displacement field. The developed 
theory is only valid for thin-walled structures. The author extended this theory to include a 
nonlinear formulation and a finite difference based numerical solution to investigate the 
structural behavior of solid orthotropic beams of arbitrary cross sections [78]. Kalfon and 
Rand [79] developed a theory for the nonlinear modeling of thin-walled composite 
helicopter blade. Inplane warping was neglected in this formulation. 
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A free vibration analysis of composite I-beams with elastic couplings under rotation 
was presented by Chandra and Chopra [80]. The developed theory was a linear analysis 
based on Vlasov theory. Constrained warping and transverse shear effects were included. 
A quasi-analytical method for the evaluation of composite box beams has also been 
developed by Smith and Chopra [59]. The cross-sectional analysis in the theory was 
performed analytically to reduce the problem to a one-dimensional beam problem which 
was solved using the finite element method. The cross-sectional analysis was based on the 
Classical Laminate Theory (CLT). The out-of-plane warping used in the formulation was 
determined using a contour analysis. Due to the many simplifying assumptions associated 
with this theory, this approach must be restricted to use as a preliminary design tool. 

Among the theories presented above, the more comprehensive anisotropic theories rely 
upon a full three-dimensional finite element solution which can become very computa- 
tionally intensive [47,48]. References [49,55-57,60-62,83,74,75,78] address 
comprehensive modeling of beams with solid cross sections. All of the closed section 
analysis procedures are based on thin-wall assumptions [50-54,58,59,63,65-73,76,77,79] 
which is not an adequate assumption when modeling the load carrying structures of tilt- 
rotor blades. It must be also be noted that when traditional beam analyses are performed in 
which the three-dimensional behavior of the structure is reduced to one-dimensional beam 
like parameters [50-54,58,59,63,65-73,76,77,79,80] (e.g. extension and three rotations), 
it is necessary to perform the cross-sectional analysis at each unique cross section. In 
many of the formulations presented above, this cross-sectional analysis may require several 
thousand degrees of freedom. As a result, this type of modeling can become quite complex 
and computationally intensive for beams which have arbitrary spanwise distributions. 
Further, many assumptions are made about the warping distributions assumed in these 
theories. As a result, the beam warping is typically determined in a manner such that the 
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average warping over the cross section is equal to zero. This, however, does not satisfy 
the condition of stress free boundary conditions on the surfaces of the beam. 

1 .4.2 Composite plate theories : In all of the previously developed closed-section beam 
models, some type of thin-walled assumption is made. However, the composite spar 
typically used in tilting prop-rotor configurations, e.g. the Advanced Technology Blades 
(ATB’s), have thicker wall sections which makes such models incapable of properly 
modeling these beams. Therefore, the need for a more general theory for adequate analysis 
of such sections is obvious. A composite box beam theory which can model sections with 
arbitrary thicknesses is, of course, more accurate as well as more realistic, since it 
eliminates all the uncertainties associated with the aforementioned assumptions. 

An objective of the present study is to develop a beam theory in which the three- 
dimensional behavior of the composite structure is not reduced to one-dimensional beam 
like parameters. By decomposing the beam into the individual walls which make up the 
beam, a beam theory is developed in which the solution approximates the exact three- 
dimensional elasticity solution. Further, the warping of the beam is determined in a manner 
such that the stress free boundary conditions on the inner and the outer surfaces of the 
beam will be automatically satisfied. To adequately model thick-walled structures, it is 
necessary to utilize an appropriate composite laminate theory in each of the walls. Several 
such theories are discussed below. 

As a first approximation to a more refined displacement field, first-order shear 
correction theories and other approximate techniques have been proposed. A first-order 
shear deformation theory (FST) was developed for general anisotropic laminated plates by 
Whitney and Pagano [81]. This theory introduces a constant shear correction term and the 
displacement field does not satisfy the necessary stress free boundary conditions on the 
upper and lower surfaces of the laminate. Despite the approximate nature of the theory, 
results for moderately thick plates were found to be more accurate than the Classical 
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Laminate Theory (CLT). Reissner [82] developed variational principles for deriving the 
plate equations from three-dimensional elasticity based on first-order shear deformation 
theories. 

A reliable six-node triangle plate/shell element was developed by Kosmatka [83] for the 
analysis of laminated composite structures based on the first-order shear deformation 
theory. Bhumbla et al. [84] used the first-order shear deformation theory to predict the free 
vibration frequencies and mode shapes of spinning laminated composite plates. The results 
from the study indicate that the Classical Laminate Theory over predicts the stiffness and 
natural frequencies. The authors later extended this work to study the buckling speed of 
spinning, laminated composite plates offset from the axis of rotation [85]. 

Rehfield and Valisetty [86,87] formulated a simple theory for the bending and 
stretching of homogeneous plates. In this theory the classical assumption that the 
transverse normal strain and the two transverse shear strains are set to zero was replaced 
with a hypothesis that the statically equivalent stresses obtained from CLT can be used to 
estimate the transverse normal strain and transverse shear strains. 

The variational asymptotical approach was used by Hodges et al. [88-90] to decompose 
the three-dimensional plate problem into a linear through-the-thickness, one-dimensional 
analysis to obtain plate elastic constants and a two-dimensional analysis to analyze plate 
deformations. Although the possibility of large deflections and rotations was considered, 
small strains were assumed in the formulation. The variational asymptotical method 
expands the strain energy in Taylor series expansion based on a small parameter, defined to 
be the thickness of the plate in this case, and as a result the application of the theory was 
confined to thin plates. 

It has long been recognized that higher-order theories are adequate for modeling 
composite laminates with arbitrary thicknesses. A third-order theory which includes the 
transverse normal stress was developed by Lo et al. [91]. The developed general theory 
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was compared to exact elasticity equations and showed excellent correlation. However, the 
retention of transverse normal shear stress terms introduces a significant amount of 
complexity to the formulation. An alternative higher-order composite laminate theory was 
developed by Reddy [92] in which the transverse shear effects were included, but the 
transverse normal stress was neglected. In this third-order theory, unlike first-order 
theories, the stress free boundary conditions were exactly satisfied on the upper and lower 
surfaces of the plate. This relatively simple theory, which has only two additional 
unknown functions from the Classical Laminate Theory or the first-order shear deformation 
theory, has been shown to be extremely accurate when compared to exact elasticity 
solutions for thick plates [93]. In a subsequent review of all third-order, two-dimensional 
theories for plates, it was shown by Reddy that all third-order theories developed over the 
last two decades are based on the same displacement field [94]. The third-order theory due 
to Reddy [92] was compared with a first-order shear deformation theory and the Classical 
Laminate Theory, for cross-ply laminates with various boundary conditions, in a study 
performed by Khdeir et al. [95]. It was found that the third-order theory out performed 
both the first -order theory and the CLT. 

Murakami {96] introduced the concept of superimposing a zig-zag linear function upon 
a FST to improve the accuracy of inplane responses. The theory was shown to be very 
accurate, but in this layerwise approach the number of unknown functions was directly 
proportional to the number of composite plies. As a result, the technique was 
computational lv expensive. This work was later extended by Toledano and Murakami [97] 
such that the piecewise linear continuous displacements were superimposed upon a 
quadratic transverse shear stress distribution. In another study performed by the authors 
[98], a zig-zag shape function was superimposed on Legendre polynomials to approximate 
the inplane displacement contributions across the plate thickness. Like the original study, 
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these theories are dependent on the number of plies and are therefore computationally 
expensive. 

A similar study was performed by Cho and Parmerter [99,100] in which a higher-order 
plate theory for composite laminates was obtained by superimposing a cubic varying 
displacement field over a zig-zag linearly varying displacement. In this theory, Heaviside 
functions were used to ensure continuity of the transverse shear stresses at the interface 
between the laminae. The results, however, were not shown to be more accurate than 
traditional third-order theories. Robbins and Reddy [101] developed a procedure for the 
modeling of thick composites using a layerwise laminate theory. The resulting layerwise 
finite element model used in this theory was capable of computing interlaminar stresses and 
other localized effects as accurately as three-dimensional finite element models. However, 
the number of degrees of freedom required in this model was approximately the same as a 
full three-dimensional finite element model and a result, the layerwise theory was 
computationally intensive. 

Chattopadhyay and Gu developed a new higher-order laminate theory for the modeling 
of delamination buckling of composite plates and shells [102,103]. Delaminations between 
layers of composite plates were modeled by jump discontinuity conditions at the 
delaminated interfaces. These discontinuities were modeled in both the lower and higher- 
order terms of displacements using Heaviside step functions. Excellent correlation with 
experimental results was obtained. An exact elasticity solution was presented by 
Chattopadhyay and Gu [104] for the buckling of simply supported orthotropic plates. The 
solutions presented in the theory indicated that the third-order theory was the most accurate 
plate theory with less than eight percent error in all cases, while the Classical Laminate 
Theory was the least accurate theory with errors of more than 100 percent for thick 


laminates. 
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1 .4.3 Extension of plate theories for beam modeling : No effort has been reported in using 
such higher-order composite plate theories in modeling composite box beams. The reason 
being the complexity of the process, especially when pre-twist and taper are included in the 
formulation as they are in the present approach. In the present research, the refined 
displacement field of Reddy [92] is used to analyze the individual walls of a composite box 
beam. As a result, the developed theory is capable of modeling box beams with arbitrary 
wall thicknesses. By decomposing the beam in this manner, the developed theory is an 
approximation to the three-dimensional elasticity solution. Therefore, there is no reduction 
of the box beam behavior to one-dimensional beam like parameters. Thus, several 
simplifying assumptions are avoided. 



2. Objectives 


The objective of this research is to develop a more general, but computationally 
efficient, theory for the adequate analysis of composite box beam sections with moderately 
thick walls. A refined higher-order displacement field is used to accurately represent the 
transverse shear stress distributions in composite laminates of arbitrary thickness which 
represent the box beam walls. The procedure developed is capable of analyzing composite 
box beam sections with pre-twist, taper and sweep to model load carrying structural 
members used in aerospace applications. Unlike previous beam theories, the present 
theory approximates the three-dimensional elasticity solution rather than reducing the beam 
properties to one-dimensional quantities. Further, the warping of the cross section in this 
theory is determined such that stress free boundary conditions are exactly satisfied on the 
inner and the outer surfaces. As a result, the model is capable of accurately describing 
thick-walled load-carrying members typically found in tilt-rotor blades such as the 
Advanced Technology Blades (ATB) on the XV-15 tiltrotor [102,103]. The developed 
beam model is general enough for applications to a wide range of wing and rotor blade 
sections. 

Next, the box beam model is used to represent the principal load carrying member in a 
tilt-rotor blade and an aeroelastic analysis is performed. The aerodynamic loading is 
obtained through an analysis based on a two-dimensional blade element momentum theory 
[104]. The coupled equations of motion are solved to determine the structural response, 
the aeroelastic stability and the aerodynamic performance of the trimmed rotor. 

Finally, a formal multidisciplinary optimization procedure is developed to investigate 
the trade-offs associated with aeroelastic stability, aerodynamic and structural design 
requirements of prop-rotors. The optimization procedure developed can be coupled with 
other analysis codes being used by the industry allowing trade-off studies to be performed 
during both preliminary and detailed design stages of the tilt-rotor aircraft development. 
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The objectives of the proposed research are summarized below. 

1 . Development and validation of a higher-order composite beam theory capable of 
modeling box beams of arbitrary shape and wall thicknesses. The theory is capable 
of modeling short aspect ratio beams with pre-twist, taper and sweep. 

2. Development of the coupled aerodynamic and structural equations of motion to 
investigate aeroelastic stability and blade structural response. 

3 . Development of a multidisciplinary optimization procedure which includes formal 
multiobjective formulation technique and an approximate analysis technique based 
on hybrid expansions. 

4. Optimization for multiple flight conditions to investigate the effects of composite 
tailoring on the overall performance of high speed prop-rotor blades. 



3. Composite Structural Modeling 


The box beam is modeled using composite laminates to represent the four walls 
(Fig. 3.1). Several different coordinate systems are used throughout this paper and are 
defined as follows. The global coordinate system (X, Y, Z) is the untwisted coordinate 
system located on the axis of rotation (Figs. 3.1 and 3.2). The global, rotated coordinate 
system (X', Y', Z') represents the coordinate system obtained by rotating the global 
coordinate system about the axis of twist by an angle 0(x). In addition, there are three 
more coordinate systems defined locally in each wall of the box beam. The local, 
untwisted coordinate system of the i-th wall is defined by (xi, yi, zj). The local, twisted 
wall coordinate system for the i-th wall which results from the global rotation (0) of the 
beam is denoted (x*, y*, z*) as depicted in Figs. 3.3 and 3.4. The introduction of sweep 

adds complexity to the formulation. For this model, sweep is defined to be in the X'-Y' 
plane (Fig. 3.5). The horizontal walls are parallel to this plane and as a result, no 



Fig. 3.1 Composite box beam. 
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additional coordinate system needs to be introduced for these walls. In the vertical walls, 
however, the sweep is normal to the walls and therefore an additional coordinate system 
must be introduced. This local, twisted and swept coordinate system is defined by (%i, r\i, 
£i). (Note that in case of the horizontal walls the coordinate system defined by (Xi, m, Ci) 
is identical to the coordinate system defined by (Xj , yj , z ; ).) Detailed explanations of the 

transformations between the coordinate systems are presented in Appendix A. 



Fig. 3.2 Beam cross section and axis of rotation. 


3.1 Refined Displacement Field 

A higher-order theory [92] is used to define the displacement field for each wall in the 
local, twisted and swept coordinate system (%, B, 0 as follows. The subscript T has 
been omitted for convenience throughout the remainder of the study. 

Qi(x> ii.O = u 0 (x.'n)+cf-^^ 1 ^ +, t , x(x. ii)l +^ 2< t | x(x-'n)+CYx(x.'n) < 31a ) 

l °x ) 

u 2 (X.il.O = Vo(X.ll) + cf-^^ + Vy(X.T|)l + C 2 <l>y(X.il) + ;Yy(X.il) (3-lb) 

v dT i J 

u 3 (X.il) = Wo(X.ll) (3 ’ lc) 
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Top wall Bottom wall 


Fig. 3.3 Horizontal wall coordinate systems. 

where u 0 . v 0 and w 0 represent the displacements at the midplane of each plate and \|/ x and 
\| f y represent the rotations of the normals to the midplane. The beam warping in each plate 
is represented by the functions <|) x , <J> y , y x and y y . The local wall deformations 
( uj , u^. u^» in the twisted coordinate system are related to local deformations (ui, U2, U3) 

in the untwisted coordinate system (x, y, z) through the following relationship. 
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or 

u = T ur ti ( 3 - 2 ) 

where T ur is the transformation matrix between the local, twisted displacements and the 
local, untwisted displacements. The local deformations in each of the plates (ui, U2, U3) 
are related to the global deformations (u, v, w) as follows (Figs. 3.1, 3.3 and 3.4). 
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Fig. 3.4 Vertical wall coordinate systems. 
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3.2 Boundary Conditions 

It is necessary to satisfy shear stress free boundary conditions at each free surface of 
the beam. This implies that in each wall, the shear stresses must be zero on both the outer 
and the inner surfaces. The imposition of these boundary conditions allow some of the 
higher-order terms to be determined in terms of some lower-order terms. This procedure, 
which is explained next, reduces the number of unknowns in the displacement field (Eqns. 
3.1). For each wall, the stress free boundary conditions are stated as 


T1, C = ±h/2) = 0, 
t), C = ±h/2 > = 0. 


(3.4) 

(3.5) 
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where h is the total wall thickness. For composite laminates made up of layers of 
orthotropic lamina, the above requirements imply that the corresponding strains must be 
zero on these surfaces. That is, 

e x ;(Z.ri,; = ±h/2) = 0, (3-6) 

e T1? (x,Ti J C = ±h/2) = 0, (3-7) 


where 


£ ^=1?- + T i = ¥x+2C<I>x+3C 2 Yx = 0, 
xs 

(3.8) 

£ ^ = dl + 9ll = Vy +2^ y + 3 C Yy - 0. 

(3.9) 

This yields the following relationships 
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3.3 Stress - Strain Relations 


Due to the fact that the stress and strain tensors are symmetric there are only six unique 


values of these quantities. Therefore, the following notation is used to define the stress and 


strain tensors in the local, untwisted coordinate system. 
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The stress and strain tensors in the local, twisted and swept coordinate system are 
expressed similarly. The generalized Hooke’s law is used to relate the stresses and the 
strains. Assuming the products of the derivatives of the displacements to be small in the 
strain formulation and using Eqn. 3.2, the strains in the local, untwisted coordinate system 


are expressed as 
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(3.14) 


In the above expressions, the short hand notations c0 and s0 are used to denote cos0 and 
sin0, respectively. It is important to note, however, that the displacement equations 
(Eqn. 3.1) are written in the local, twisted and swept coordinate system. The relationship 
between the derivatives in the local, untwisted coordinate system and the local, twisted and 


swept coordinate system are written as 
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— = cos 0^ -sin 0^7, (3.15b) 

3y oTl oL, 

A = sine^ + cose47. ( 3 - 15c ) 

8z oTj OL, 


Details of the transformation relationship between the coordinate systems are found in 
Appendix A. Using Eqns 3.1, 3.14 and 3.15, the strains are written as follows. 
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(3. 16e) 


(3.16f) 


{ _W 0,„ + Vy - ^-C 2 ¥y JcosG Z * 0 > 


From the above equations, the following relationship between the local strains in the 
untwisted and the twisted coordinate systems is obtained. 


£ i - T ij( £ j + 9,xPj z o, x ^j 


(3.17) 


where 8° is the strain in the local, twisted coordinate system in the absence of pre-twist 
and sweep, jlj is the additional strain components due to pre-twist and 9 x is the twist rate. 
The additional strain components due to the sweep are denoted and the sweep rate is 
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z* It must be noted that additional sweep terms appear only in the vertical walls due to 

°,x 

the definition of sweep used in the modeling (see Fig. 3.5). The total strain in the local, 
untwisted coordinate system (x, y, z) is denoted 8j and Ty is the transformation matrix 
between the strains in the local, twisted coordinate system and the strains in the local. 


untwisted coordinate system. This transformation matrix is expressed as follows. 


T- 


'1 0 0 0 0 0 ‘ 

0 c 2 0 s 2 0 — c0s0 0 0 

0 s 2 0 c 2 0 c0s0 0 0 

0 2c0s0 -2c0s0 c 2 0 - s 2 0 0 0 

0 0 0 0 c0 s0 

0 0 0 0 -s0 c0 


(3.18) 


Note that the transformation matrix for the strains ( Ty) is different from the transformation 


matrix for the displacements (T ur ). The local inplane strains in the absence of pre-twist or 


sweep are derived as follows, 
e? = e? + £k{ + C 3 k? 

S2 = e 2 + C k 2 + ? K l 

£6 = e 6 + + CM 


(3.19a-c) 


The out of plane strains are expressed similarly as 


£3=0, 

S^eS+CM. ( 3 ' 19d - f) 

E§=e§ + C 2 K5- 

The additional strains due to beam pre-twist are as follows 

P-! = n? + ?m- 1 + CV 2 + cV? + A A • 

M-2 = fe = P-4 = 

P -5 = ^5 + 

ji 6 = H6 + CM6 + C 2 |4 + + cVt 


(3.20a-f) 
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The additional strain components due to the sweep are 

= ^ 3 = 0 , 0 4 =0, 0 5 = 0, (3.21a-f) 

d 6 = dg+C 2 «6. 

where the nonzero components of the individual strains are described in Appendix B. 


3.4 Energy Formulation 

The beam equations of motion are derived using Hamilton’s principle [93] which 

assumes the following form. 
l 2 

8 j*(U -T + W e )dt = 0 (3.22) 

l l 


where §( ) represents the variation of ( ) and U, T and W e represent the total beam strain 
energy, kinetic energy and external work, respectively. Using variational principles, Eqn. 
3.22 may be rewritten in terms of the individual plate quantities as follows. 


l 2 f N 

J5> 


-STi+SWe. 


\ 

dt = 0 


(3.23) 


q v i=l 2 

where N is the total number of walls (N = 4 for a box beam). The individual strain energy 
density (U 0 ) in each plate is calculated as follows. 

£ i 



(3.24) 


0 


Using the generalized Hooke’s Law (a t = Qij£j), Eqn 3.24 is rewritten as 

U 0 = ^-QijEiej. 0-25) 


where repeated indices (i, j = 1, 2, — , 6) indicate summation and cq is the strain tensor. 
The quantities Qij denote the full three-dimensional material properties in the local. 
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untwisted coordinate system which, for laminates made of orthotropic plies, is represented 
as follows. 


Qii = 


Qn 

Ql2 

Ql3 

0 

0 

Ql6 

Ql2 

Q22 

Q23 

0 

0 

Q26 

Ql3 

Q23 

Q33 

0 

0 

Q36 

0 

0 

0 

Q44 

Q45 

0 

0 

0 

0 

Q45 

Q55 

0 

_Ql6 

Q26 

Q36 

0 

0 

Q66_ 


(3.26) 


The material properties in terms of the local, untwisted coordinate system (x, y, z) are 
written in terms of the material properties in local, twisted and swept coordinate system 
(X, rt, 0 as follows. 

O • = T 0 T • (3.27) 

where Q mn represents the material properties in the local, twisted coordinate system. The 
total strain energy in the i-th wall (Uj) is then written, using Eqns. 3.18 and 3.26 as 
follows. 


i=Ju 0i dV 

V 

— ~ J (^m + ®,xA-m — z o x ^m)Tik^miQmn^njTjl^n + ® xM-n — z o ?x ^n 


(3.28) 


* ~ 


V 


where V indicates integration over the volume of the wall. Due to the orthogonality of the 
transformation matrix, T , (T ik f mi = 6^, where 5km is the Kronecker delta) Eqn. 3.28 

is simplified as follows. 

h/2 

Ui = — J + ^xM-m ~ z O )X ^m)Qmn^n + ^xP-n “ z o >x ^m) (3.29) 

Q- h/2 


where dO is the differential area (d£2 = dxdr)). The strain energy can be rewritten using 
Eqns. 3.19-3.21 and 3.29. 
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Ui 


iJ 


m 


£! 


-A-mn ®mn ^mn ^mn ^mn 

®mn ^ran E mn F mn G mn 

^mn ^mn ^mn ^mn ^mn 


G 


H, 


0 T 


^mn A mn VJ mn “mn w mn 

^mn G mn H mn O mn P mn j 


B n da 


(3.30) 




Q. 

where m, n = 1, 2, •••, 6 and the laminate stiffness matrices (A - P) are defined in each of 
the walls as follows. 


h/2 


(A,B,D,E,F,G,H, 0 ,P) = J Q(l,« 2 ,^^^^ 7 ,£ 8 )df; 


(3.31) 


-h/2 


For composite laminates, Q is ply dependent and as a result it is a function of the thickness 
coordinate Therefore, this matrix cannot be taken out of the integral in Eqn. 3.31. The 
vector B m is defined as 




( e m + Fm^,x — z o x ^m) ( K m + Fm^,x) ( K m + Fm^,x z o x ^m) 
( K m“^Fm^,x) Pm^,x] 


(3.32) 


The external work due to applied loads and body forces (W e ) in the i-th wall is written 
as 

W e . = J* XjUjdV + J tjUjd 5 j = 1 , 2 , 3 (3.33) 

V 5 

T 

where Uj is the displacement vector defined as [u^ U 2 U 3 ] , Xj, X 2 , and X 3 are the 

body forces in the X, Y and Z directions, respectively. Applied surface tractions over the 
region of the surface 5 are denoted t h t 2 and t 3 , along the respective directions. 

The total kinetic energy in the i-th wall is expressed as 

T i = “ J p Vj VjdV j = 1., 2, 3 

v 


(3.34) 
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where Vi is the velocity vector defined as 
3u; 

v j = — ~ + Qk x iq , 

3t 


(3.35) 


Uj is the displacement vector, O is the rotational velocity about the Z axis and r\ is a 
position vector from the axis of rotation to an arbitrary point in the i-th wall (Fig. 3.6). The 


position vector is written as follows. 


r = X 0 i + Y 0 j + (% + u 1 )e x + (ti + tj 0 + u 2 )q^ + (C + C 0 + u 3 )e^, (3.36) 

where X 0 and Y 0 are the offsets in the global, untwisted coordinate system from the axis of 
rotation to the center of twist and r| 0 and are the local offsets to an arbitrary point on the 
wall expressed in terms of the local, twisted and swept axis system. Using Eqns. 3.35 and 
3.36 and the coordinate transformations defined in Appendix A, the kinetic energy is 
expressed as follows. 

<| u [ 2 + 2uj ^ + Co + w)Qti “ 2G i + Ti* + u 2 )n^ + (C + Co + u 3 ) Q* 
v 

- 2(C + C* +u 3 )(ri + Tlo + S 2)flri Q ^ + ( T l + T lo + Q 2) 

. ^ 2 1 (3.37) 

- ■!u; 2 + 2 u 2 (x + Xo + Q l) Q ^ + (x + Xo + 5 l) a ^| 

+ ;D ? 2 -2a 3 (x + Xo+ a l) a Ti + (x + Xo + Q l) a ^| dv 


where the rotational and position vectors (Eqns. 3.35 and 3.36) are rewritten as 
r = (x + X<> + £i i )c y + (ti + T|o + “2 )e n + (C + Co + «3 f; > 

Ok = Q^e^ + fTe-, 

and (') represents a differentiation with respect to time. 


(3.38) 

(3.39) 
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Fig. 3.6 Rotational offsets and position vector. 


3.5 Incremental Stiffness Due to Rotation 

To properly obtain the natural frequencies of rotating plates and beams, it is necessary 
to consider the stiffening effects which arise due to the stresses generated by rotation. This 
procedure involves two steps. First the stresses of the beam are calculated for centrifugal 
forces only. Then the incremental stiffness due to rotation is computed and is added to the 
original stiffness matrix. The procedure for calculating the incremental stiffness matrix is 


outlined below. 

The additional strain energy due to centrifugal stiffening in each of the walls (U c f) is 
given by [105] 

u cfl 4Jb 2 ^ + < a '/ + ^ 2 (°* + ^) 


v 


(3.40) 


-200^(0^^ - 20^0^0^ - 2co % co dV 


where a> x , and oo £ are the rotations about %, rj, £ axes, respectively. It is important to 
note that in the above equation, the strains are not included in the formulation based on the 
assumption that their effects are small compared to the rotations. Equation 3.40 may be 
written in matrix notation as follows. 
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u cfi4JK “’l ®d 


V 


0 TJ 

a xn 

_a xC 

®x 


a x 





( a % + °ri)_ 

CO^ 


dV 


(3.41) 


where the stresses Gij are due purely to the centrifugal force in each individual wall. The 
rotations are calculated as follows. 


co v = — 

K o 


a>n = - 
11 2 


co r = — 


9u 3 

3u 2 ^ 


a; J 

f d“i 

9u 3 ^ 

Us 

dX J 

f 952 

3ui ") 


(3.42a-c) 


dx 

Using Eqns. 3.12 and 3.42 the rotations are rewritten as 


co 




1 drj 


0W O 

>T l “ 

3X 

1 

f 3v o 

'5 "2 

[dX 


V h 
1 (. 4 ) 


1 9 

h 2 ) 


V y> 


3u 0 ^ 


(3.43a-c) 


3.6 Variational Method 

The variation of strain energy is written as follows. 

5Uj = J(i^Q. mn 5B n d£2 m, n = 1, 2, — , 6 (3.44) 

n 


where the variation of the strain vector B m is expressed as the sum of the variation of the 
strain in the untwisted coordinate system, the variation of the strain due to pre-twist and the 
variation of the strain associated with sweep as follows. 



Slim = [(8e° + 5|i° e * - 8< z ; x ) (64 + 5nLe, x ) 

(84+54e x -8*m z 0,x) ( 8K m +5 Hm e ,x) ^l 0 ,x] 
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(3.45) 

The variation of the potential energy of the applied loads is expressed as 

8W e . = J XjSujdV + J tjSSjdS, j = 1, 2, 3. (3.46) 

V 5 

Due the presence of the rotational velocity, the variation of the kinetic energy can be 
expressed in terms of four individual components as follows. 

STi = 5T m , + 6T ci + 5T ki + 5T fi (3.47) 

where the subscript ‘m’ denotes the component of the kinetic energy that is used to form 
the mass matrix. Similarly, the subscripts ‘c\ ‘k’ and ‘f indicate the components that are 
used in the formulation of the gyroscopic matrix, the stiffness matrix and the forcing 
vector, respectively. The exact form of these variations are as follows. 

5T m . = J* P [ U]8uj + U28U2 + U38U3I dV (3.48) 

v 

8T Cj = Jp | (2u : £2~ - 2630^)88 ! - 2C]Q^8u 2 + 2u,£2 ri 8u 3 dV (3.49) 

V 

8 T k . = Jp [( + Q = ju 1 5u 1 + (ogu 2 - Qti^u 3 )5u2 

v (3.50) 

- ( -f2 r ,Q-u 2 +n^u 3 j5u 3 J dV 

5T f . = Jp[ (Q“ +^)(x + Xo)Su 1 + (^(ri + Ti;)-aT 1 Q c C)Su2 

v (3.51) 

+ (-n^(ri + rio) + a^o)5u3 dV 
The quantity 8T mi , in each wall, is rewritten as 
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A4 
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i 

Ajk 

c OT 1 
<1 

§4 


Ajk 

Ajk 

Ajk ' 

§4 ■ 

Q 

A] k 

Ajk 

i 

54 


j, k=l, 2, 3 (3.52) 


where X°, X 1 and A, 3 correspond to the zeroth, the first and the third-order components 
(in 0 of the displacement field, respectively. These quantities are defined as follows. 

T 

A.°=[u 0 v 0 w 0 J 



(3.53a-c) 


The density matrices, Aj k , are defined as 

h/2 

(A^ k ,A 1 jk ,A] k ,A] k ,Aj k ,A^ k )= Jp(U,C 2 ,C 3 .C 4 .C 6 )dC. (3-54) 

-h/2 

Although simple closed form expressions cannot be obtained for the remaining components 
of the variation of the kinetic energy, they can be computed in a straight forward manner as 
discussed later. 

The variation of the strain energy due to the centrifugal stiffening in each wall is written 
as 

SU cfj = [< X mn 5co n dV, m, n = 1, 2, 3 (3.55) 

V 


where 
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®m > 


dw ° | 1 

dx 2 


f 


4 rl 


^ F* 


if 9v o 9u 0^ 


(3.56) 


and 


t 


mn 


®r\ 

1 

Q 

3 


_a %ri 

a x 


.“°x5 


(°x +a 


(3-57) 


It must be noted that the stresses in Eqn. 3.57 are due to centrifugal forces only which are 
determined based on the forces associated with Eqn. 3.51. Therefore, Eqn 3.51 is used 
only to determine the steady state stresses due to rotation. After these stresses are 
calculated and the incremental stiffnesses are determined, the forcing terms associated with 
this equation are no longer included in the formulation. 


3.7 Solution Procedure 

The solution of the equations of motion is obtained using a two-dimensional finite 
element formulation in the local, twisted and swept coordinate system of each individual 
plate (%, T), 0- A four noded plate element is used to discretize the individual plates of the 
beam. This element is C 1 continuous in the zeroth order displacements (u G , v 0 , w Q ) and is 
C° continuous in the higher order terms (\j/ x , \|/ y ). As a result, the element contains 11 
degrees of freedom per node which are defined in terms of the nodal degree of freedom 
vector as follows. 


q = 


3u n <3u n 3v n 5v 0 5w 0 dw 0 


-iT 


’ 3ti ’ °’ 3x ’ an ’ °’ 9x ’ 9 ti 


(3.58) 
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3.7.1 Continuity conditions : To maintain the continuity of displacements throughout the 
entire beam, constraints are imposed at the comers of each individual plate as follows. 

1 u 0 (x>'n = bi) = 2 u 0 (x. , n = o) 

1 v 0 (x>'n = b l) = - 2 w 0 (X>'n = °) (3.59a-c) 

1 w 0 (x. , n = b i) = 2 v o(x.'n = °) 

where the preceding superscripts ‘1’ and ‘2’ denote walls 1 and 2, respectively and bi is 
the width of wall 1 (Figs. 3.3 and 3.4). It must be noted that these equalities must be 
satisfied for all values of %. Therefore, the partial derivatives of the above equalities, with 
respect to %, must also be satisfied. To ensure that the angle between the walls remains 90° 
after deformation, the following constraints are imposed on the rotations about the X-axis. 

1 w 0>tl (X^ = b l) = 2 w 0j11 (X,T 1 = 0) (360) 

V y (x.'n = b i) = 2 Vy(x.-n = o) 

Similar sets of constraints are derived at each of the four comers of the box beam. 


3.7.2 Finite element formulation : The finite element approach is used to solve the 
complete beam equations of motion (Eqn. 3.23). Denoting q as the nodal degree of 
freedom vector, it is possible to express the strain (Eqn. 3.5) in the following form 
£ i =(r 1J +e, x o iJ -z: x H 1J )q J . (3.6i) 

The partial derivatives of the strain vector with respect to qj are then written as 


de.j 

a qj 


~ r ij + G.x^ij z Ov^ij- 


(3.62) 


Note that quantities T, <F and E can be expanded in terms of £ as follows 

r^ = \r§ rj rl i*].{i C C 2 C 3 C 4 } 

= .{i c C 2 C 3 C 4 }. 

-[*§ *|j C C 2 C 3 C 4 } 


(3.64) 
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and 


“ij 


=N ^H 1 ?2 } 

= ^*0 c 2 }- 


(3.65) 


The rotation terms (coi) may also be expressed as a function of the nodal degrees of 
freedom as follows 


C0i - 0jjqj 


(3.66) 


so that 
OcOj 


a qj 


= 0y. 


Similarly, the displacement vector u may be written as 
s, 


U 1 

G 2 

[U 3 


’ u j 

s vj ^ 


’w 


JJ 


(3.67) 


(3.68) 


or 




(3.69) 


such that partial derivatives of the displacements with respect to qj are as follows 

^I = Sii . (3.70) 

3qj 


Note that the displacement matrix can similarly be written in terms of the thickness 


coordinate as 


Sij=N s!j sfj].{l q C 3 } 
= *v{i C C 3 }- 


(3-71) 


Using the finite element approach, the coupled (dynamic/aerodynamic) beam equations 


of motion are written in matrix form as follows 
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M q(t) + C q(t) + K q(t) = F(t) (3.72) 

where M, C and K are the mass, the damping and the stiffness matrices, respectively and 
F is the external force vector representing loads corresponding to the beam nodal degrees 
of freedom (q). Note that in Eqn. 3.72, the C matrix represents the gyroscopic (Coriolis) 
effects and not damping of the system. As a result, the equations of motion as formulated 
represent an undamped system. Therefore a proportional damping (two percent) is 
assumed in the model to represent the structural damping. The proportional damping is 
determined based on the natural frequencies obtained through consideration of the mass and 
the stiffness terms only in Eqn. 3.72. These terms are augmented to the C matrix. It must 
be noted that additional terms analogous to damping and stiffness terms will be introduced 
through aerodynamic loading. These terms are augmented to the appropriates matrices in 
Eqn. 3.72 as explained in detail in Chapter 5. The matrices and forcing vector due to the 
structural contribution only are expressed as follows. 



and the forcing vector is 



(3.76) 
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N 

f cf = XI 
i=l 


Jp(q^(s u {%+Xo} + S wCoj ^o| s v^o 

V 

+ ^(s u {x + Xo} + s v{t| + %}) 


(3.77) 


dV. 


Note that the matrix C in Eqn 3.75 does not include the proportional damping which can be 
determined only after an eigenanalysis of Eqn. 3.72 is performed. The matrix denoted K* 
(Eqn. 3.75) corresponds to the stiffness matrix obtained before the addition of the 
incremental stiffness due to rotation. Using Eqns 3.75 and 3.76 to determine the stresses 
due to beam rotation, the incremental stiffness matrix is calculated as follows. 


N 


K 


cf 


=x 


i=l 


0 T i 0 dV 


(3.78) 


where the stress matrix (X ) is determined from the solution to 


K* q = f c f ( 3 - 79 ) 

The total stiffness matrix (K) used in Eqn 3.72 is now written as a combination of K* and 

K c f as follows 


K = K* + K cf . 


(3.80) 



4. Composite Beam Results and Validation 


The mesh sizes necessary for both individual plate and complete beam analyses are 
determined by performing a detailed convergence study. The details of this convergence 
study are presented in Appendix C. Next, beam correlations are presented in order to 
demonstrate the adequacy of the individual wall elements. Finally, the beam model (both 
thin-walled and thick-walled sections) is correlated with available experimental results and 
results obtained using a general purpose finite element code. Details of the correlation 
study are presented in the following sections. In the following sections, the elemental 
mesh size is defined as M x N where M is the number of chordwise elements and N is the 
number of span wise elements. For the beam model, a consistent mesh is used in every 
wall and the mesh size presented corresponds to an individual wall. A finite element 
representation of a beam with a 10 x 4 mesh is illustrated in Fig. 4.1. 



4.1 Higher-Order Plate Verification Studies 

To demonstrate the importance of including the transverse shear terms in the 
formulation for the individual plates and to prove how well the present theory can capture 
these effects, results are first presented for individual plates. The accuracy of the plate 
theory is established from these validation studies, which are presented below. 
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4.1.1 Square, untwisted plate : Results are presented for a simply supported square, 
orthotropic plate under uniform loading (Fig. 4.2). The material properties of the plate are 
listed in Table 4.1. Figures 4.3 and 4.4 present the variations of the normalized center 
deflection, the normal axial stress (Gi) and the transverse shear stress (G5) with plate 
thickness. The results of the present theory are compared with the those obtained using an 
exact elasticity approach [106]. Note that in Fig 4.4, results using the classical laminate 
theory (CLT) are not presented because in case of the axial stress, they are nearly identical 
to the results from the higher-order theory and in case of the transverse shear stress, they 
are zero. The figures indicate that the higher-order plate theory correlates very well with 
the exact elasticity solution [106]. Also, for moderately thick to very thick plates, the 
assumption of zero transverse shear stresses in CLT can introduce significant errors 
(Fig. 4.4). 



Fig. 4.2 Orthotropic square plate with fixed boundary conditions. 
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Table 4. 1 Summary of orthotropic material properties 105 


Ei = 20.83 x 10 6 p.s.i., E 2 = 1.094 x 10 6 p.s.i., 
G12 = 6.10 x 10 6 p.s.i., G13 = 3.71 x 10^ p.s.i., 
p 12 = 0.44 



Width to thickness ratio, a/h 

Fig. 4.3 Normalized center deflection of fixed orthotropic square plate under uniform 

distributed load. 



Fig. 4.4 Normalized stresses of fixed orthotropic square plate under uniform distributed 

load. 
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4.1.2 Plate with twist : From the convergence study it is determined that a 4 x 12 mesh 
(665 degrees of freedom) provides converged results and therefore this mesh size is used in 
all analyses involving the twisted plates. To validate the accuracy of the higher-order 
theory, the results are compared with published results. The first four natural frequencies 
for the plate are presented in Table 4.2. The natural frequencies are nondimensionalized as 
follows. 


= CDj 



(4.1) 


D = Eh 3 /l2(l-v 2 ). ( 4 ’2) 

The NASTRAN [105] results presented are due to MacBain [107] which uses a mesh size 
of 1 1 x 24 to yield a total of 1265 degrees of freedom. The Rayleigh-Ritz solution is due 
to Barton [108] in which an 18 term expansion for the deflection is used. In the table and 
following figures ‘F’ is used to denote a flexural or bending mode, ‘T’ is used to denote a 
torsional mode and ‘PM*’ is used to indicate the i-th plate mode. From the table it is noted 
that although both the NASTRAN and the present approach correlate very well with the 
Ritz solution, the present approach is more accurate and requires fewer total degrees of 
freedom. 


Table 4.2 Nondimensional frequencies of a flat, cantilevered plate 
(v = 0.3, L/w = 2.33) 


Mode 

Ritz 

NASTRAN 

% Difference* 

Present 

% Difference* 

IF 

3.47 

3.43 

- 1.2 

3.43 

- 1.1 

IT 

17.10 

16.74 

- 2.1 

16.87 

- 1.3 

2F 

21.58 

21.36 

- 1.0 

21.45 

- 0.6 

2T 

55.00 

53.66 

- 2.4 

54.08 

- 1.7 


* Percent difference from Ritz solution. 
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The natural frequencies of the twisted, isotropic plates over a range of tip twist values, 
are presented in Figs. 4.5 - 4.7. In addition to the results obtained using the present 
approach, experimental data [107] and results obtained using NASTRAN are also 
presented. In these plates, the twist is assumed to vary linearly along the span. Results are 
presented for tip twist values of 0 t = 0°, 12°, 17°, 23.5°, 30° and 38° to correlate with 
available experimental data. In addition, results using the present approach are also 
calculated at pre-twist values of 45° and 60° to further examine the trends. From the 
figures, excellent correlation is noted in all cases, with a possible exception of the first plate 
mode (PMi) for low tip twist values. However, it must be noted that the experimental 
value for this particular mode is also not available [107]. Using the present approach, its 
value is very close to the third torsional mode (3T, Fig. 4.7). The proximity of these 
natural frequencies might explain why it was possible to experimentally determine only one 
of these values (3T). This also suggests that the current theory does accurately predict the 
natural frequency of this mode and the results obtained using NASTRAN are under 
predicted. Overall, the present approach yields more accurate results than NASTRAN. 
This is due to the fact the elements used in NASTRAN are based on first-order shear 
deformation theories which only approximate the transverse stress. 
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Fig. 4.7 Natural frequency as a function of tip twist (modes 7 - 10). 

4. 1 .3 Vibration of thick, swept plates : To validate the higher-order theory for thick plates 
with sweep, the first eight nondimensionalized frequencies are calculated for several 
different plate configurations. The plate geometry is shown in Fig. 4.8 and plates with two 
different length-to-width ratios corresponding to a thick plate (a/b = 5) and a very thick 
plate (a/b = 2) are investigated. In both cases, three different width-to-thickness ratios 
(b/h = 0.5, 1 and 2) are used. The results obtained using the present approach are 
compared with two different numerical results which were presented by McGee and Leissa 
[109]. The first set of results in Ref. 109 were obtained using a three-dimensional Ritz 
solution in which a 6 x 4 x 4 (288 degrees of freedom) mesh was used to discretize the 
displacement field. In addition, results were also obtained using NASTRAN solid 
elements (CHEXA) for a mesh size of 14 x 14 x 3 (2520 degrees of freedom). It must be 
noted that both of these results represent truly converged results. In the present approach, 
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a mesh size of only 6 x 4 (269 degrees of freedom) is used which also represents 
converged results. Graphical results for the thick plate (a/b - 5) with a unit width-to- 
thickness ratio (b/h = 1) are presented in Figs. 4.9 - 4.16. The results from the other five 
cases are presented in Tables. 4.3 - 4.7. 



Fig. 4.8 Definitions of swept plate. 


The first natural frequency of the thick plate with unit width-to-height ratio is shown in 
Fig. 4.9 from which several important observations are made. It is seen that increasing the 
sweep angle will increase the natural frequency. Also, all three techniques are in excellent 
agreement. The second natural frequency is presented in Fig. 4.10 and indicates that for 
zero sweep angle all three approaches are in very good agreement. With increases in the 
sweep angle, the natural frequency in case of the Ritz solution increases more rapidly than 
the other two techniques. This due to the fact the stress free boundary conditions on the 
edges of the plate are exactly satisfied in the Ritz solution. In both the NASTRAN results 
and the present approach these boundary conditions are not satisfied. In cases of moderate 
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sweep (up to 30°), the difference between the present approach and the Ritz solution is less 
than three percent. However, despite the fact that the NASTRAN solution involves an 
order of magnitude increase in the total degrees of freedom, the present approach is a better 
approximation to the Ritz solution. 

Similar trends are observed in the higher-order modes (Figs. 4. 1 1 - 4.16). In particular 
it observed that in general the natural frequency increases more dramatically with sweep in 
the Ritz solution than with the either the NASTRAN results or the present approach. In 
some cases, the natural frequency is slightly over predicted in the present approach for zero 
sweep. However, in all cases as the sweep angle increases the results using the present 
approach are as good or are better than NASTRAN despite the fact that the number of 
degrees of freedom used in NASTRAN is an order of magnitude larger than that used in the 
present approach. In addition to the results presented in Figs. 4.9 - 4.16, the results 
presented in Tables 4.3 - 4.7, for various plate thicknesses and length-to-width ratios, also 
show similar trends. From all of these results, the adequacy of the higher-order theory to 
model very thick, swept plates is demonstrated. 
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Fig. 4.9 Comparison of the natural frequencies of the first mode for thick, swept plates. 



Fig. 4.1 1 Comparison of the natural frequencies of the third mode for thick, swept plates. 
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Fig. 4.10 Comparison of the natural frequencies of the second mode for thick, swept 

plates. 
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Fig. 4.13 Comparison of the natural frequencies of the fifth mode for thick, swept plates 
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Fig. 4.14 Comparison of the natural frequencies of the sixth mode for thick, swept plates. 


o 

£ 

c/) 



j0;0uiBJBd i(0U0nb0Jij 


Fig. 4.15 Comparison of the 








56 


Table 4.3 Natural frequency parameters of thick, swept cantilevered plates; 

a/b = 0.5, b/h = 5 


Analysis 

CP 

Sweep Angle 
15° 30° 

45° 

Mode 1 

Ritz 

3.1227 

3.2166 

3.5501 

4.1425 

Present 

3.1054 

3.2161 

3.5459 

4.0818 

NASTRAN 

3.1238 

3.2217 

3.4885 

3.8582 

Mode 2 

Ritz 

4.2261 

4.3542 

5.0753 

7.1414 

Present 

4.3052 

4.4372 

4.9293 

6.2007 

NASTRAN 

4.2821 

4.4073 

4.9183 

6.2958 

Mode 3 

Ritz 

6.8552 

7.2038 

8.2730 

10.2007 

Present 

6.8642 

6.9203 

6.9706 

6.5989 

NASTRAN 

6.7974 

6.8581 

6.9365 

6.6924 

Mode 4 

Ritz 

8.0642 

8.3654 

9.3562 

11.4351 

Present 

7.5768 

7.7075 

8.1924 

9.4582 

NASTRAN 

7.3213 

7.4760 

8.0192 

9.2653 

Mode 5 

Ritz 

12.6494 

12.5656 

13.3447 

16.6312 

Present 

12.5927 

12.1570 

11.6274 

11.7331 

NASTRAN 

12.5705 

12.0459 

11.4246 

11.5358 

Mode 6 

Ritz 

12.9342 

12.9723 

14.7501 

18.9970 

Present 

13.0389 

13.1159 

13.5085 

15.1500 

NASTRAN 

12.9816 

12.3409 

12.7782 

13.3430 

Mode 7 

Ritz 

13.5239 

14.3488 

17.1083 

22.5717 

Present 

13.3894 

13.6536 

14.4809 

15.7681 

NASTRAN 

13.3094 

13.1870 

13.5041 

14.2523 

Mode 8 

Ritz 

13.9620 

14.7390 

17.7805 

24.5030 

Present 

13.9822 

14.1885 

15.2488 

16.4272 

NASTRAN 

14.4315 

14.3078 

15.1933 

16.7757 


Ritz solution 6x4x4 mesh; present solution -6x4 mesh; MSC/NASTRAN CHEXA 
values - 14 x 14 x 3 mesh. 
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Table 4.4 Natural frequency parameters of thick, swept cantilevered plates; 


a/b = 2, b/h = 5 


Analysis 

OP 

Sweep Angle 
15° 30° 

45° 

Mode 1 

Ritz 

3.3397 

3.3432 

3.0183 

2.9961 

Present 

3.4077 

3.4750 

3.6810 

4.0514 

NASTRAN 

3.4114 

3.5002 

3.7805 

4.3144 

Mode 2 

Ritz 

12.4593 

12.5382 

12.7095 

. 11.6690 

Present 

13.7012 

13.9136 

13.8973 

12.8721 

NASTRAN 

13.2833 

13.5085 

13.8810 

12.9673 

Mode 3 

Ritz 

14.3907 

14.4765 

15.4606 

16.0299 

Present 

14.4754 

14.3496 

14.6975 

16.4236 

NASTRAN 

14.4521 

14.3258 

14.4257 

16.6212 

Mode 4 

Ritz 

19.5960 

19.7466 

20.6890 

23.6116 

Present 

20.4807 

21.2286 

23.5735 

28.1333 

NASTRAN 

20.3647 

21.1291 

23.4390 

27.4966 

Mode 5 

Ritz 

38.7711 

37.5107 

35.3326 

35.659 

Present 

43.3142 

42.6169 

42.2367 

44.2161 

NASTRAN 

41.9086 

41.6010 

42.1593 

46.0781 

Mode 6 

Ritz 

52.2825 

49.9046 

51.9786 

64.1992 

Present 

52.2598 

51.1111 

49.4186 

47.6519 

NASTRAN 

52.3346 

50.6122 

48.6333 

47.2957 

Mode 7 . 

Ritz 

52.4692 

50.3584 

55.3036 

68.8087 

Present 

54.4937 

57.2897 

59.9279 

61.6180 

NASTRAN 

53.4834 

55.9763 

59.9545 

62.1115 

Mode 8 

Ritz 

54.8140 

54.0790 

57.9550 

71.5920 

Present 

55.8318 

57.3844 

64.8321 

75.3250 

NASTRAN 

54.7376 

56.9011 

62.7143 

72.3215 


Ritz solution 6x4x4 mesh; present solution -6x4 mesh; MSC/NASTRAN CHEXA 
values - 14 x 14 x 3 mesh. 


Table 4.5 Natural frequency parameters of very thick, swept cantilevered plates; 

a/b= l,b/h = 2 
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Sweep Angle 


Analysis 

OP 

15° 

30° 

45° 

Mode 1 
Ritz 
Present 
NASTRAN 

2.9463 

2.9402 

2.9397 

3.0096 

3.0113 

3.0017 

3.3164 

3.2317 

3.1809 

3.8412 

3.6310 

3.4554 

Mode 2 
Ritz 
Present 
NASTRAN 

4.4178 

4.3996 

4.3957 

4.6232 

4.4303 

4.4244 

4.6768 

4.4928 

4.4802 

6.5555 

4.4518 

4.4421 

Mode 3 
Ritz 
Present 
NASTRAN 

5.1815 

5.4081 

5.1470 

5.3635 

5.5332 

5.2527 

6.1756 

5.9641 

5.6351 

7.6914 

6.9015 

6.4876 

Mode 4 
Ritz 
Present 
NASTRAN 

10.5391 

10.4674 

10.5200 

9.9399 

10.3013 

10.2640 

11.3566 

9.9520 

9.8067 

13.3667 

9.5993 

9.4300 

Mode 5 
Ritz 
Present 
NASTRAN 

10.9792 

11.4304 

10.7864 

10.645 

11.6097 

10.8493 

13.5265 

12.044 

10.9455 

18.6951 

12.3342 

10.8145 

Mode 6 
Ritz 
Present 
NASTRAN 

11.7535 

11.9244 

11.6626 

11.3436 

12.1834 

12.0043 

13.5265 

12.7421 

12.6506 

18.6951 

13.0822 

13.0694 

Mode 7 
Ritz 
Present 
NASTRAN 

14.4674 

15.0781 

14.3273 

12.5292 

14.4062 

13.7176 

14.3425 

14.3705 

13.6732 

18.9918 

15.3392 

14.0382 

Mode 8 
Ritz 
Present 
NASTRAN 

16.1660 

15.1308 

14.4794 

14.4180 

16.3263 

15.3277 

15.9530 

17.8745 

16.5283 

25.3940 

18.7010 

16.7953 


Ritz solution 6x4x4 mesh; present solution -6x4 mesh; MSC/NASTRAN CHEXA 
values - 14 x 14 x 3 mesh. 
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Table 4.6 Natural frequency parameters of very thick, swept cantilevered plates; 

a/b = 0.5, b/h = 2 




Sweep Angle 


Analysis 

CP 

15° 

30° 

45° 

Mode 1 
Ritz 
Present 
NASTRAN 

2.2304 

2.2483 

2.2375 

2.2987 

2.3041 

2.2757 

2.5790 

2.4744 

2.3791 

3.1205 

2.6396 

2.5398 

Mode 2 
Ritz 
Present 
NASTRAN 

2.7039 

2.7457 

2.6884 

2.8662 

2.7681 

2.7535 

3.3293 

2.7883 

2.7861 

4.1110 

2.7899 

2.6920 

Mode 3 
Ritz 
Present 
NASTRAN 

2.7577 

2.7908 

2.7289 

2.8973 

2.8774 

2.7675 

3.4345 

3.1797 

3.0452 

4.8256 

3.8379 

3.4531 

Mode 4 
Ritz 
Present 
NASTRAN 

4.7438 

4.6172 

4.2550 

4.8997 

4.6900 

4.2921 

5.3477 

4.6509 

4.3260 

6.6697 

4.6932 

4.3008 

Mode 5 
Ritz 
Present 
NASTRAN 

5.0722 

5.0371 

5.0096 

5.0376 

4.8628 

4.8003 

5.4416 

4.9193 

4.5566 

6.8077 

5.1200 

4.5895 

Mode 6 
Ritz 
Present 
NASTRAN 

5.5855 

5.5929 

5.5375 

5.8947 

5.6647 

5.1834 

6.6160 

5.5440 

5.0861 

8.1640 

5.9195 

5.3187 

Mode 7 
Ritz 
Present 
NASTRAN 

5.5932 

5.6144 

5.5698 

5.9298 

5.6754 

5.8741 

7.1205 

6.0995 

5.9334 

9.7608 

6.5709 

5.4364 

Mode 8 
Ritz 
Present 
NASTRAN 

5.8219 

5.7794 

5.7046 

6.1000 

5.9215 

5.7498 

7.3500 

6.4047 

6.1511 

10.0930 

6.9770 

6.5378 


Ritz solution 6x4x4 mesh; present solution -6x4 mesh; MSC/NASTRAN CHEXA 
values -14x14x3 mesh. 



Table 4.7 Natural frequency parameters of very thick, swept cantilevered plates; 

a/b = 2, b/h = 2 
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Sweep Angle 


Analysis 

CP 

15° 

30° 

45° 

Mode 1 
Ritz 
Present 
NASTRAN 

3.2545 

3.2684 

3.2667 

3.2823 

3.3259 

3.3220 

3.2943 

3.5032 

3.4865 

3.3134 

3.8267 

3.7563 

Mode 2 
Ritz 
Present 
NASTRAN 

5.7869 

5.7902 

5.7970 

5.8978 

5.7398 

5.7466 

5.9928 

5.5589 

5.5692 

6.6404 

5.1488 

5.2041 

Mode 3 
Ritz 
Present 
NASTRAN 

9.9090 

10.5107 

9.9270 

9.8656 

10.6933 

10.0899 

10.4693 

11.3024 

10.6543 

11.0825 

12.5047 

11.8026 

Mode 4 
Ritz 
Present 
NASTRAN 

16.3311 

16.8113 

16.4537 

16.3082 

17.1731 

16.7024 

17.4151 

18.2583 

17.3983 

19.7762 

19.0607 

18.3928 

Mode 5 
Ritz 
Present 
NASTRAN 

20.9459 

20.9039 

20.9648 

20.9449 

20.4445 

20.2626 

21.8063 

19.7675 

19.4645 

25.7577 

20.0736 

18.9263 

Mode 6 
Ritz 
Present 
NASTRAN 

21.9266 

22.3327 

21.8993 

23.563 

22.9537 

22.7818 

25.7634 

23.9712 

24.0112 

32.3171 

24.6472 

24.8581 

Mode 7 
Ritz 
Present 
NASTRAN 

29.5524 

31.4764 

29.5814 

28.8472 

31.2740 

29.6120 

30.5295 

31.3982 

30.1386 

24.1152 

33.0256 

31.9981 

Mode 8 
Ritz 
Present 
NASTRAN 

37.2420 

39.2395 

37.4149 

30.8320 

40.5095 

38.1270 

34.8060 

43.6541 

39.7678 

45.4330 

40.7105 

41.1186 


Ritz solution 6x4x4 mesh; present solution -6x4 mesh; MSC/NASTRAN CHEXA 
values - 14 x 14 x 3 mesh. 
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4.2 Static Thin- Walled Beam Verification Study 

To validate the developed procedure, correlations are made with available experimental 
results on a thin-walled box beam [1 10] and a previously developed analytical model [59]. 
The analytical model developed in Ref. 59 is a one-dimensional thin-walled beam model, 
based on the Classical Laminate Theory (CLT). In this approach, the out-of-plane warping 
effects are determined using on a contour analysis. The details of the beams studied are 
presented in Table 4.8. The cross ply and symmetric beams are all subjected to two 
different loading conditions, a 1 lb. bending load at the tip and a 1 lb.-in. tip moment. The 
anti-symmetric beams are subjected to a 1 lb. axial load at the tip and a 1 lb.-in. tip 
moment. 

4.2.1 Cross plv : The bending slope of the cross ply beam under a 1 lb. tip bending load is 
presented in Fig. 4.17, which compares the experimental data [110], the results of the 
quasi-analytical model [59] and the results from the present study. Further, results from a 
beam finite element model reported in Ref. 59 are also presented. As mentioned in Ref. 
59, this two-dimensional finite element technique is as accurate as a full three-dimensional 
finite element model for the particular beams studied in that report. From Fig. 4. 17, it is 
seen that all three modeling techniques under predict the bending slope of the cross ply 
beam. This can possibly be attributed to errors in the experimental results arising from 
fiber alignment problems that are typically encountered during fabrication and curing of the 
beam. A slight shift from the desired (0° and 90°) ply orientations can introduce additional 
coupling terms which will reduce bending stiffnesses, thereby increasing the bending 
slope. Overall, there is still good correlation between all three modeling techniques in this 
case. The twist angle of the cross ply beam due to a 1 lb.-in. tip moment is shown in Fig. 
4.18. From the figure, excellent correlation is observed for both the quasi-analytical model 
[59] and the present approach. To further investigate the possible coupling with slightly 
shifted angle plies, the bending slope and elastic twist for the same loading conditions are 
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calculated for a “cross ply” [6784 °] 3 beam. From the figures it is seen that both the 
bending slope (Fig. 4.19) and the elastic twist (Fig. 4.20) now correlate extremely well 
with the experimental results. 

Table 4.8 Details of composite beams 59 ’ 1 10 


Flanges Webs 

Top Bottom Left Right 


Cross Ply 

[0790°] 3 

[0790°] 3 

[0790°] 3 

[0°/90°] 3 

Symmetric 15° 

[15°] 6 

[15°] 6 

[15%15°] 3 

[157-15°] 3 

Symmetric 30° 

[30 o ] 6 

[30°] 6 

[307-30°] 3 

[30%30°] 3 

Symmetric 45°- 

[45 °] 6 

[45 °] 6 

[45%15°] 3 

[45°/-45 °] 3 

Anti-symmetric 15° 

[15°] 6 

[-15 c k 

[15°]6 

H5°]6 

Anti-symmetric 30° 

[0730°j3 

[0%30°] 3 

[0°/30°] 3 

[0%30°] 3 

Anti-symmetric 45° 

[0745°] 3 

[0%45 °] 3 

[0°/45 °] 3 

[07-45°] 3 

Length = 30 in., width = 

= 0.953 in., depth 

= 0.53 in., ply thickness 

= 0.005 in, number of plies = 6, wall 

thickness = 0.030 in. 

Mechanical properties: E L = 20.59 x 

10° p.s.i., E t = 

1.42 x 10° p.s.i.. 

G lt = 0.89 x 10 6 p.s.i.. 

v LT = 0.42. (Cross ply dimensions: width 

= 2.06 in., depth = 

1.025 in.) 


A Experimental 



Fig. 4.17 Bending slope of cross ply beam under 1 lb. tip bending load. 



0.0005 


A Experimental 



Fig. 4.18 Twist angle of cross ply beam under 1 Ib.-in. tip moment. 



Fig. 4.19 Bending slope of “cross ply” [6°/84°]3 beam under 1 lb. tip bending load. 



Fig. 4.20 Twist angle of “cross ply” [6°/84°]3 beam under 1 lb.-in. tip moment. 
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4.2.2 Symmetric beams : The bending slope of the symmetric 15° beam under a 1 lb. tip 
bending load is presented in Fig. 4.21 where good correlation is observed between both 
modeling techniques. There are two sets of experimental data presented in this figure (as 
well as in the following figures) due to the fact that two separate beams were tested in Ref. 
1 10. The induced twist due to tip bending load is presented in Fig. 4.22, which shows that 
the present approach slightly over predicts the twist angle at the tip compared to the quasi- 
analytical model of Ref. 59. Overall, both models show good correlation with the 
experimental data. 



Fig. 4.21 Bending slope of symmetric 15° beam under 1 lb. tip bending load. 



Fig. 4.22 Bending induced twist of symmetric 15° beam under 1 lb. tip bending load. 
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The bending slope of the symmetric 30° beam is presented in Fig. 4.23. This figure 
also shows the results from the beam finite element method presented in Ref. 59. From the 
figure it is observed that the quasi-analytical technique slightly under predicts and the 
present approach slightly over predicts the bending slope. The beam finite element method 
[59] over predicts the slope more significantly than the current approach. Once again, good 
correlation exists between all techniques. The induced twist due to bending load is 
presented in Fig. 4.24. The quasi-analytical technique again under predicts the response. 
The results from the present study correlate extremely well with the experimental data in 
this case. 


A Experimental 



Fig. 4.23 Bending slope of 30° symmetric beam under 1 lb. tip bending load. 



0.020 n 
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Experimental 



Fig. 4.24 Bending induced twist of 30° symmetric beam under 1 lb. tip bending load. 


The results from the symmetric 45° beam subjected to a 1 lb. tip bending load are 
presented in Figs. 4.25 and 4.26. As shown in Fig. 4.25, the bending slope is slightly 
over predicted by the present model although the correlation with experimental data is still 
very good. The quasi-analytical technique [59] also correlates well. In case of the induced 
twist due to the bending load, the trends are significantly different (Fig. 4.26). The quasi- 
analytical method greatly under predicts the twist in this case while the present approach 
correlates extremely well with the experimental data. Further, the results are in excellent 
agreement with those predicted by the variational asymptotical approach (VABS) due to 
Cesnik et al. [69]. The approach reduces the cross-sectional properties into one- 
dimensional beam properties based on an expansion in terms of a small parameter which is 
defined to be the beam height divided by the beam length. The theory also includes both 
inplane and out-of-plane warping and is well suited for thin-walled box beams. 
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Fig. 4.25 Bending slope of [45°]6 thin-walled beam under 1 lb. bending load at tip. 


A Experimental 



Fig. 4.26 Bending induced twist of [45°]6 thin-walled beam under 1 lb. bending load at 

tip. 

The results from the 1 lb.-in. tip moment loading case, for all three symmetric beams, 
are presented in Figs. 4.27 and 4.28. Since the variation of the response is linear, only the 
results at the mid span location (x/R = 0.5) are presented. From Fig. 4.27 it is seen that 
there is scattered correlation with the experimental bending slope, but it must be noted that 
the actual values of the slope are very small (on the order of 0.0005 rad). In case of the 15° 
symmetric beam, all three approaches over predict the slope when compared to the 
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experimental values. In case of the symmetric 30° beam, the present approach and the 
beam finite element method again over predict the slope, whereas the quasi-analytical 
method [59] correlates well with the experimental data. Finally, in case of the symmetric 
45° beam, the quasi-analytical method greatly under predicts the slope whereas the beam 
finite element method and the present approach correlate very well. Overall, however, all 
three techniques show good correlation with the experimental data, particularly when the 
actual magnitude of the slope is taken into consideration. The comparisons of the twist 
angles are presented in Fig. 4.28. In case of the symmetric 15° beam, all three techniques 
slightly over predict the twist angle. However, in case of the symmetric 30° beam and the 
symmetric 45° beam, both the present approach and the beam finite element method 
correlate extremely well with the experimental data while the quasi-analytical technique 
significantly under predicts this behavior. 


□ Experimental 
Quasi-Analytical 


Beam F.E.M. 
Present 



Fig. 4.27 Twist at x/R = 0.5 for 1 lb.-in. tip moment of symmetric beams. 



I I Experimental 
^ Quasi-Analytical 

ESI Beam F.E.M. 
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Fig. 4.28 Bending slope at x/R = 0.5 for 1 lb-in. tip moment of symmetric beams. 

By examining Figs. 4.17 - 4.28, it is observed that the quasi-analytical model of Ref. 
59 correlates well with experimental data for the cross ply beam and the 15° symmetric 
beam. In case of the 30° symmetric beam, the quasi-analytical model begins to under 
predict the behavior. This is particularly evident in Fig. 4.22 which presents the induced 
twist due to a 1 lb. tip bending load. The quasi-analytical method greatly under predicts the 
behavior of the symmetric 45° beam in all cases studied with the exception of the bending 
slope due to a 1 lb. tip bending load. It is also observed from these figures that although 
the present approach slightly over predicts the beam behavior for lower ply angles, the 
technique correlates very well with the experimental data in cases with larger ply 
orientations. Further, in cases where the present approach does over predict the behavior 
of the beam, compared to the experimental values, the beam finite element method [59] 
exhibits similar trends (Figs. 4.17, 4.27 and 4.28). 
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4.2.3 Anti-symmetric beams : The results from the three anti-symmetric beams (Table 4.8) 
are presented in Figs. 4.29 and 4.30. Since the response is linear, only the results at the 
mid span location are presented. As in Figs. 4.27 and 4.28, the actual magnitude of the 
twist angle, must be noted in these figures. The twist angle due to a 1 lb.-in. tip moment is 
presented in Fig. 4.29. In case of the 15° anti-symmetric beam, the quasi-analytical model, 
the beam finite element method and the present approach all correlate very well the 
experimental data. For the other two anti-symmetric beams, all three approaches predict 
very similar results, slightly under predicting the response compared to the experimental 
results. Similar trends are observed in Fig. 4.30 which presents the twist angle for a 1 lb. 
axial tip load. For the anti-symmetric 15° beam all three approaches slightly over predict 
the twist. However, all three approaches predict the behavior very well for the anti- 
symmetric 30° beam and the anti-symmetric 45° beam. 

Overall, both the quasi-analytical model of Ref. 59 and the present approach correlate 
well with the experimental data in case of the anti-symmetric beams. Both techniques 
predict the same trends for all three beams. The present approach offers small 
improvements over the quasi-analytical results for the anti-symmetric 15° beam. 
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□ Experimental 
Quasi-Analytical 


H Beam F.E.M. 



Fig. 4.29 Twist at x/R = 0.5 for 1 ib.-in. tip moment of anti-symmetric beams. 


I I Experimental 
f^| Quasi-Analytical 
[5*3 Beam F.E.M. 



Fig. 4.30 Bending slope at x/R = 0.5 for 1 lb. tip axial load of symmetric beams. 
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4.3 Static Thick-Walled Beams Results 

To demonstrate the importance of including transverse shear effects in the beam 
formulation, results are now presented for a series of thick-walled beams. Due to the lack 
of available experimental data, only numerical results are presented. Two different 
composite lay-ups are used which correspond to the symmetric 15° beam and the symmetric 
45° beam (Table 4.8). The beams studied have a length-to-width ratio (L/c) of 2.5 and a 
width-to-height ratio (c/d) of 2. Since the goal is to investigate the effects of thick-walled 
beams, the wall thicknesses used in this study are 0.25 in. resulting in the values of the 
width-to-thickness ratio, c/h = 8, in the horizontal walls and the height-to-thickness ratio, 
d/h = 4, in the vertical walls. These two beam configurations are subjected to a 100 lb. 
bending load at the tip as well as a 100 lb.-in. tip moment. 

Figure 4.31 presents the elastic twist for the thick-walled 15° symmetric beam subjected 
to a tip bending load. From the figure, it is observed that in addition to the fact that the 
local twist in the four individual walls is nonlinear, the average twist in each of the four 
walls is also nonlinear. This is different from the trends observed in case of the thin-walled 
beams where the local twists in the individual walls are nearly identical. In general, the 
local values of twist are very close to the average values which are presented in Fig. 4.27. 
Thus the average twist is a good representation of the beam twist for thin-walled beams. In 
the thick-walled case, however, the local twist differs significantly in the individual walls 
and as a results it is difficult to designate a value of “beam twist” for the entire cross 
section. Figure 4.31 also shows that the twist in the vertical walls are not equal to each 
other, while the twist in the horizontal walls are equal. This is due to the fact that the ply 
angles in the horizontal walls are all +15° (measured in the global, twisted coordinate 
system) and the two opposite walls are therefore mirror images of each other in their local 
coordinate systems. In the vertical walls, however, the stacking sequence in the opposite 
walls differs by a sign change (e.g. ±15° in the right wall and +15° in the left wall). This 
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nonlinear trend is more observable in Figs. 4.32 and 4.33 which presents the induced twist 
of the symmetric 45° beam subjected to a tip bending load. In addition to the results 
obtained using the present approach, results obtained using NASTRAN are also presented. 
The NASTRAN results are obtained using QUAD4 plate elements which are elements 
based on first-order shear deformation theory. From Figs. 4.32 and 4.33, good correlation 
is observed between the present approach and NASTRAN in each of the walls. The 
nonlinear behavior of the vertical walls is more dominant in this case. The average twist 
for both approaches is presented in Fig. 4.34 where the correlation is again noted. In 
general, NASTRAN slightly over predicts the results compared to the present approach. 
This is due to the fact that only first-order shear deformation effects are included in the 
NASTRAN elements (QUAD4) used in this study. 



Fig. 4.31 Bending induced twist of thick-walled 15° beam under 100 lb. tip bending load. 



Twist angle xlO 3 (rad) oq‘ Twist angle xlO 3 (rad) 
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4.32 Bending induced twist of thick-walled 45° beam under 100 lb. tip bending load; 

horizontal walls. 



Right (Present) 

0 Right (NASTRAN) 

Left (Present) 

--G- Left (NASTRAN) 


Fig. 4.33 Bending induced twist of thick-walled 45° beam under 100 lb. tip bending load; 

vertical walls. 
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Fig. 4.34 Comparison of average bending induced twist of thick-walled 45° beam under 

100 lb. tip bending load 

A more complete explanation why the twist varies in each of the individual walls is 
presented in Figs. 4.35 - 4.37. Figure 4.35 shows the application of the tip bending load 
to the individual walls of the beam. In Fig. 4.36, a schematic diagram of the resulting 
displacements for unconnected walls (that is, individual plates) is shown where it is seen 
that displacements in the horizontal walls are described by a translation and a rotation. This 
is due to the fact that the plies in these walls are all of the same value (e.g. +45°) and 
therefore the laminate is unbalanced. Since both the vertical walls comprise balanced 
laminates (e.g. ±45°), there is no rotation in these walls and the displacement is purely 
translational. From Fig. 4.36 it is seen that the rotational displacement of the horizontal 
walls restricts the translational motion of the left vertical wall whereas the rotational motion 
in the horizontal walls is complimentary to the translational motion of the right wall. As a 
result, the horizontal walls become cambered which in turn induces large shearing stresses 
in the vertical walls. This type of bending behavior is shown in Fig. 4.37 which presents 
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the midplane deformation (greatly magnified) of the tip cross section for the 45° symmetric 
beam subjected to a tip bending load obtained using both the present approach and 
NASTRAN. Both techniques show shearing of the vertical walls and the cambering of the 
horizontal walls. Careful examination of Fig. 4.37 also reveals that the walls do not remain 
perpendicular to each other (at the comers) after deformation in case of the NASTRAN 
results. This is due to the fact that no such constraints are imposed in the NASTRAN 
formulation. In the present formulation, constraints are imposed on the rotations at the 
comers and the walls do remain perpendicular after deformation. This explains the increase 
in camber of the horizontal and the decrease in the transverse shear of the vertical walls 
when compared to the present results. 



Fig. 4.35 Schematic of load distribution under tip bending load. 



Fig. 4.36 Schematic of individual wall displacements under tip bending load. 
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Fig. 4.37 Deformation of cross section at tip due to tip bending load; 45° beam. 


The twist angle due to a 100 lb.-in. tip moment for the 15° symmetric, thick-walled 
beam is presented in Fig. 4.38. Similar trends as those obtained with the tip bending load 
are observed. In this case, however, the average twist of all four individual walls is very 
nearly linear. In comparison, in case of the thin-walled, symmetric 15° beam, the twist 
angle due to a 1 lb.-in. tip moment is linear in each of the four walls, except at the tip where 
the values are slightly larger for the horizontal walls compared to the vertical walls. The 
trends for the thick-walled beams are more observable in Figs. 4.39 and 4.40 which 
presents the twist angles of the thick-walled 45° symmetric beam due to a 100 lb.-in. tip 
moment for both the present approach. These figures also present the corresponding 
NASTRAN results. A comparison of the average twist using both techniques is presented 
in Fig. 4.41. Good correlation is observed from these figures. The small differences are 
once again attributed to the fact that QUAD4 elements used in NASTRAN only includes 
first-order shear deformation effects and the lack of constraints at the comer rotations in the 
NASTRAN formulation. 

The results presented above indicate the importance of including transverse shear 
effects in the beam formulation. Further, the results also show that in general the twist of a 
composite beam is a local quantity which can be defined only locally, at a point. The 
definition of the twist at the centroid of a beam (or some other arbitrary point) is an 
approximation which in case of thick-walled sections can be erroneous. 
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Fig. 4.38 Twist of thick-walled 15° beam under 100 lb.-in. tip moment. 



Fig. 4.39 Elastic twist of thick-walled 45° beam under 100 lb.-in. tip moment; horizontal 

walls. 
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Spanwise coordinate (in) 


Fig. 4.40 Elastic twist of thick- walled 45° beam under 100 Ib.-in. tip moment; vertical 

walls. 



Spanwise coordinate (in) 


Fig. 4.41 Comparison of average elastic twist of thick-walled 45° beam under 100 lb.-in. 

tip moment. 
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4.4 Dynamic Results 

In following sections, dynamic results for several different beam configurations are 
presented, including both thin-walled and thick-walled beams. Due to the lack of 
experimental data available, the results obtained using the present method are compared 
with those obtained using NASTRAN. In case of the composite beams, QUAD4 elements 
(first-order shear deformable) are used. In case of the thick-walled isotropic beam, both 
CHEXA (solid elements) and QUAD4 elements are used. 


4.4.1 Thin-walled beams : The natural frequencies of the first several modes are calculated 
using the present approach and NASTRAN. These frequencies which are 
nondimensionalized as follows 


Xj = C0j 


IphLT 

D 


(4.3) 


D = Ej h 3 /l2(l-v 12 2 ). (4-4) 

are presented in Table 4.9. The beam bending and chordwise bending modes are denoted 
‘B’ and ‘C\ respectively, the torsional modes are denoted £ T’ and extension modes are 
denoted ‘E\ An element mesh size of 4 x 30 (each plate) is used in the present approach 
and since the NASTRAN elements are linear, a mesh size of 12 x 30 is used in each plate. 
From Table 4.9, very good correlation is observed between the predicted frequencies 
obtained using the present approach and NASTRAN. In general, the natural frequencies 
are slightly higher in case of the present model. This is once again due to the fact that a 
third-order displacement field is used in the present approach and the NASTRAN elements 
use only a first-order displacement field. Further, fairly significant differences observed in 
the torsional frequencies, with the present model being more stiff torsionally. The higher 
natural frequencies predicted by the present approach (compared NASTRAN results) is due 
to the constraints imposed on the comer rotations which are not present in the NASTRAN 
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formulation. It must also be noted that the present approach predicts a purely inplane 
warping mode which occurs before the fifth chordwise bending mode and the second 
torsional mode. This mode is not captured with NASTRAN. 


Table 4.9 Natural frequencies of symmetric, thin-walled 45° beam 



Frequency Parameter 

Modes 

Present 

NASTRAN 

(4x10 mesh) 

(QUAD4, 12x30 mesh) 

B1 

27.57 

27.05 

Cl 

47.79 

46.98 

B2 

172.32 

169.12 

C2 

297.23 

292.48 

B3 

480.13 

472.39 

C3 

822.92 

812.08 

B4 

930.75 

921.98 

T1 

1458.03 

1301.71 

B5 

1498.28 

1514.79 

C4 

1581.34 

1570.51 

El 

1735.09 

1718.58 

Warp 

2363.50 

N.A. 

C5 

2553.36 

2568.33 

T2 

3354.10 

3388.66 


To demonstrate the effect of inplane and out-of-plane warping on beam dynamic 
deformation, several mode shapes are presented for the 45° symmetric composite beam 
studied in Ref. 110. Due to the stacking sequence of this composite beam (Table 4.8), 
flap-lag coupling is absent. However, both bending-torsion coupling and extension-shear 
coupling are present. Further, since the beam consists of only ±45° plies, it is extremely 
rigid in torsion and the first torsional natural frequency (C0ti) is 52.8 times the fundamental 
natural frequency (0) o = 18.83 Hz). As a result, there is no warping, either inplane or 
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out-of-plane, in the first five beam bending modes (C0b5 = 54.3 coj) and in the first four 
chordwise bending modes (0) c4 .= 57.4 co 0 ). To illustrate the lack of coupling and/or 
warping in the first several modes, the fourth chordwise bending mode is shown in 
Fig. 4.42. (In these figures the dots denote the original undeformed shape of the beam.) 
However, there exists a purely inplane warping mode which occurs before the second 
torsion mode ((% = 122 ©o)- This mode, whose natural frequency is 85.7 times the 
fundamental frequency, is illustrated in Fig. 4.43. It must be noted that this warping mode 
is not predicted by NASTRAN. The lack of significant out-of-plane warping for this 
composite beam is due to its thin-walled construction. Due to the very thin walls, the beam 
bending motion is accounted for by pure bending in the horizontal walls and the chordwise 
bending motion is accounted for by pure bending in the vertical walls. The bending 
motions are much larger than any inplane shear in the walls and as a result there is very 
little warping. 


3-D view 
Z 



Fig. 4.42 Fourth chordwise bending mode of [45°]6 thin-walled beam. 
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3-D view 
Z 



top view 


Fig. 4.43 Inplane warping mode of [45°]6 thin-walled beam. 

4.4.2 Pre-twisted thin-walled beams : In addition to the beams studied in Refs. 59 and 
110, the frequencies and mode shapes are also presented for the symmetric, 45° beam with 
a 30° twist from root to tip. Table 4.10 presents the nondimensionalized frequency 
parameters obtained using both the present approach and NASTRAN. From the table it is 
seen that there is good correlation between the two approaches, although the present 
approach predicts slightly higher values in general. Similar trends as those of the untwisted 
beam are noted in this case. The differences are once again attributed to the different 
elements being used in the formulation. The current approach models the transverse strains 
more accurately than NASTRAN and exactly satisfies the stress free boundary conditions 
on the inner and outer surfaces of the beam. As is the case for the untwisted beam, a 
purely inplane warping mode is predicted using the present approach. This mode which 
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occurs before the second torsion mode (and fourth chordwise mode) is not predicted by 
NASTRAN. 

Interesting mode shapes for the pre-twisted beam are presented in Figs. 4.44 - 4.46. 
Unlike the untwisted beam which displays very little coupling between bending modes, the 
pre-twisted beam starts to exhibit coupling between beam bending and chordwise bending 
as early as the first two modes (CDbi = 0) o and co cl = 1.56 co 0 , a> 0 = 20.5 Hz). This 
coupling is more significant in the second bending modes which are presented in Fig. 4.44 
(co b2 = 5.79 co 0 ) and Fig. 4.45 (co C 2 = 9.54 ©o) despite the fact that their natural 
frequencies are not close to each other. Similar to the untwisted beam, a purely inplane 
warping mode exists for the pre-twisted beam with a natural frequency of 77.4 times larger 
than the fundamental frequency as shown in Fig. 4.46. This mode which occurs before the 
second torsional frequency (co t2 = 109 co 0 ) is again not captured by NASTRAN. 
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Table 4.10 Natural frequencies of symmetric, thin-walled 45° beam with 30° pre-twist 



Frequency Parameter 

Modes 

Present 

NASTRAN 

(4x10 mesh) 

(QUAD4, 12x30 mesh) 

B1 

30.60 

27.10 

Cl 

47.66 

46.64 

B2 

177.17 

170.72 

C2 

292.07 

287.89 

B3 

485.96 

479.92 

C3 

811.72 

788.35 

B4 

941.58 

948.29 

T1 

1410.57 

1295.00 

B5 

1483.31 

1457.52 

C4 

1595.27 

1627.26 

El 

1735.09 

1719.61 

Warp 

2367.11 

N.A. 

C5 

2731.88 

2598.25 

T2 

3354.22 

3341.20 
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3-D view 
Z 



top view 


Fig. 4.44 Second beam bending mode of [45°] 6 thin-walled, pre-twisted beam. 


3-D view 
Z 



top view 


Fig. 4.45 Second chordwise bending mode of [45°] 6 thin-walled, pre-twisted beam. 
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3-D view 
Z 

Y 

side view 


EKEDIlZim 

: top view 

Fig. 4.46 Inplane warping mode of [45°] 6 thin-walled, pre-twisted beam. 

4.4.3 Thick-walled beams : In addition to the beams studied in Refs. 59 and 110, natural 
frequencies and mode shapes are also presented for a thicker and shorter version of the 
beam with two different sets of material properties. Complete details of the thick-w ailed 
beams studied are listed in Table 4.11. The first set of material properties corresponds to 
an isotropic beam. The second beam is made of orthotropic laminae with identical lay-up 
and material properties as the symmetric 45° beam listed in Table 4.8. 

Table 4.1 1 Details of moderately thick beam 

Length = 10 in., width = 2 in., depth = 1 in., 
ply thickness = 0.0667 in., number of plies = 6, 
total wall thickness = 0.4 in. 

Isotropic material properties 
E = 10 x 10 6 p.s.i., v = 0.3, p = 0.1 lb m /in 3 
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The natural frequency parameters for both the present approach and for NASTRAN 
results are presented in Table 4.12. Since this beam has thicker walls, solid elements 
(CHEXA) are used to model the beam in NASTRAN. A 12 x 30 x 3 mesh is used for a 
total of 17280 degrees of freedom. For comparison purposes, QUAD4 elements (12 x 30 
mesh, 4320 degrees of freedom) are also used to calculate the frequencies. A 4 x 10 
element mesh is used (1896 degrees of freedom) in the present approach. Good correlation 
is again observed between the present approach and NASTRAN CHEXA results although 
far fewer degrees of freedom are used in the present approach. The NASTRAN QUAD4 
elements also correlate well with both the present approach and the solid (CHEXA) 
elements. In general, the solid elements predict slightly higher values of the bending 
frequencies and slightly lower values of the torsional and extensional frequencies when 
compared to the other two techniques. The present approach predicts slightly higher values 
compared to the QUAD4 elements. The differences in the torsional frequencies are again 
caused by the lack of rotation constraints at the corners of the beam in both NASTRAN 
results. For this beam, as was the case of thin-walled beams, the present approach does 
predict warping modes which are not predicted by either of the NASTRAN formulations. 
However, these modes occur only after the fifth torsional mode (T5) and as a result they 
are not presented in Table 4.12. 

In the isotropic beam case, it is seen from Fig. 4.47 that out-of-plane warping is 
present as early as in the second chordwise bending mode (co c2 = 10.5 C0 o , 0) o = 313 Hz). 
This is due to the presence of inplane shear in the side walls. The third torsional mode 
(co t3 = 34.8co 0 ) is presented in Fig. 4.48. A careful examination of this figure shows a 
small amount of both inplane and out-of-plane warping. This warping, which represents 
somewhat of a three-dimensional camber, is greatest near the node points. The camber 
effect is due to the shearing of the cross section. 



Table 4.12 Natural frequencies of thick-walled isotropic beam 


Frequency Parameter 


Modes 

Present 
(4x10 mesh) 

NASTRAN 
(CHEXA, 
12x30x3 mesh) 

NASTRAN 
(QUAD4, 
12x30 mesh) 

B1 

8.27 

8.57 

8.24 

Cl 

16.30 

16.53 

16.24 

B2 

46.41 

49.28 

46.28 

T1 

74.29 

60.75 

61.18 

C2 

86.73 

88.21 

86.23 

B3 

114.74 

124.02 

114.34 

El 

130.22 

126.56 

130.04 

•pn 

197.03 

179.82 

178.06 

B4 

205.08 

208.83 

196.04 

C3 

216.61 

215.51 

203.33 

T3 

288.08 

292.11 

279.82 

B5 

289.18 

316.45 

286.29 

C'4 

299.52 

345.00 

334.98 

e: 

356.81 

378.49 

367.09 

C5 

432.30 

488.40 

472.97 
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To investigate the effect of wall thickness on the warping of composite beams, natural 
frequencies and mode shapes are calculated for a thicker and shorter composite beam with 
cross-sectional dimensions that are approximately twice that of the symmetric 45° beam 
originally studied in Ref. 1 10. The length of this beam is one-third the previous beam. 
The laminate stacking sequence and material properties are identical to the previous beam 
(see Tables 4.8 and 4.11). 

The natural frequency comparisons between the present approach and NASTRAN are 
presented in Table. 4.13 and show good correlation between the two techniques. It must 
be noted that although this beam does represent a thick-walled beam, QUAD4 (plate) 
elements are used. This is due to the fact that if solid elements (CHEXA) were used, each 
unique lamina would require a separate element. For this simple beam, that would require 
a discretization into 6 elements through the thickness of each wall. For more complex 
beams, this number would increase making the approach computationally expensive. Also 
numerical conditioning is an issue since very thin solid element do not behave very well 
[105]. Similar trends as those obtained in the previous three beams are again observed in 
this beam. In particular, the present approach predicts slightly higher natural frequencies 
for the lower modes and slightly lower natural frequencies for the higher modes compared 
to those predicted using NASTRAN. This is again due to the differences in the 
displacement field and the constraints on the comer rotations. 

As in the previous beam, for the first several modes there is no coupling between the 
beam bending, the chordwise bending and the torsional modes. However, for this thicker 
beam, a slight out-of-plane warping effect is observed in the first chordwise bending mode 
(co cl ), whose natural frequency is only 2.10 times larger than the fundamental frequency 
(C0 o = 224 Hz). The second chordwise bending mode (co C 2 = 13 - 5 ©o) is shown in Fig. 
4.49. A significant amount of out-of-plane warping is observed in this mode 
predominantly due to the shearing in the upper and the lower walls. Figure 4.49 also 
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shows that this mode is uncoupled, although its natural frequency is close to the first 
torsional mode (co tl = 11.5 (D 0 ). 


Table 4. 13 Natural frequencies of thick-walled symmetric 45° beam 


Frequency Parameter 


Modes 

Present 
(4x10 mesh) 

NASTRAN 
(QUAD4, 12x30 mesh) 

B1 

2.77 

2.70 

Cl 

5.82 

5.72 

B2 

16.67 

16.28 

T1 

31.76 

30.25 

C2 

37.32 

31.25 

B3 

43.85 

43.08 

El 

44.16 

43.47 

B4 

78.42 

77.27 

C3 

79.76 

77.48 

T2 

96.48 

82.70 

B5 

102.50 

114.11 

T3 

115.14 

117.63 

E2 

117.29 

127.22 

C4 

125.29 

132.64 

T4 

128.73 

145.54 
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Fig. 4.49 Second chordwise bending mode of [45°]6 thick-walled beam. 


There is, however, a significant amount of coupling between the fourth beam bending 
mode (C 0 b 4 = 28-3 0) o ) and the third chordwise bending mode (co C 3 = 28 - 8 ®o)- To 
illustrate this coupling, the third chordwise bending mode is shown in Fig. 4.50. This 
coupling, which is clearly due to the fact that their natural frequencies are very close, 
causes a significant amount of both inplane and out-of-plane warping as depicted in Fig. 
4.50. Unlike the second chordwise bending mode which remains uncoupled from the 
nearby first torsion mode, the third chordwise bending mode is slightly coupled with the 
second torsion mode (co t2 = 34.8 co 0 ) as well as with the fourth beam bending mode. 
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3-D view 



side view 



top view 


Fig. 4.50 Third chordwise bending mode of [45°]6 thick-walled beam. 


Finally, to illustrate the importance of including both inplane and out-of-plane warping 
in the beam formulation, the second extensional mode (co e2 = 42.3 co 0 ) is presented in Fig. 
4.51. A significant amount of warping in observed in this mode which is slightly coupled 
with the third torsion mode (co t3 = 41.6 (0 o ) and is largely coupled with the fourth 
chordwise bending mode (G ) C 4 = 45.2 co 0 )- Unlike the previous modes which have 
primarily linear warping, both the inplane and the out-of-plane warping in this mode are 
nonlinear. Of particular interest is the “necking” effect observed near the cantilevered edge 
shown the side view and the nonlinear out-of-plane cross-sectional camber which is 
demonstrated in the top view. This three-dimensional warping is a result of the shearing 
effects and are significant due to the thick-walled construction of the beam. 




95 


3-D view 
Z 



top view 


Fig. 4.5 1 Second extensional mode of [45°] 6 thick-walled beam. 



5. Aerodynamic Modeling 


The aerodynamic formulation is based on the two-dimensional compressible 
aerodynamic representation developed by Smith [1 1 1] and later modified by Talbot [1 12] 
for the formulation of axial flow performance analysis. The modifications by Talbot 
include an empirical correction to the two-dimensional stall behavior to represent the high 
lift capability demonstrated by rotors and propellers. In the initial study performed by 
Smith [1 1 1] an empirical fit was performed on NACA 63 and 64 series airfoil families in 
order to supply a functional relationship between maximum lift coefficient and sectional 
thickness and camber for incompressible flow. Detailed expressions for the coefficients of 
lift, drag and pitching moment (ci, Cd and c m ), which represent the high lift capability of 
rotary wings in post stall angle of attack region, are found in Ref. 106. These functional 
relationships were later modified by Talbot [112] to model the Advanced Tiltrotor Blades 
(ATB) [113,114]. These relationships were then modified by McCarthy et al. [13] to 
include blade sweep. A similar algorithm was proposed to model the post stall delay due to 
rotation by Corrigan and Schillings [115]. In this study, the formulation is extended to 
include blade dynamic effects. The coefficients of lift, drag and the pitching moment, 
obtained using this analysis, are presented in Figs. 5.1 - 5.3 for a typical section of the 
Advanced Technology Blade over a range of Mach numbers. 


5.1 Aerodynamic Loads 

The blade element theory used in the algorithm is due to Glauert [116]. In this 
formulation, the sectional lift and drag are resolved into elemental thrust and torque for each 
section of the blade. The force and momentum equations for thrust and torque assume the 
following form. 

dTj = 47trp(V 00 + Ui )t)jdr , (5.1) 


dT 2 = — pW 2 c(q cos A - c d sin A)dr , 


(5.2) 
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dQ] = 4ra* 2 p(V 00 +Dj)u T dr, 


dQ 2 = ^ pW 2 c(cj sin A + c d cos A)rdr , 


(5.3) 

(5.4) 


where dT, dQ and dr represent the section thrust, torque and element length, respectively, 
Voc is the forward velocity, and u T represent the inflow and swirl velocities, respectively 
and W is the magnitude of the resultant velocity. The chord length and radial locations are 
denoted c and r, respectively and p is the air density. The quantities c\ and c d represent the 
sectional coefficients of lift and drag, respectively. 



Fig. 5.1 Sectional lift coefficient distribution; Advanced Technology Blade. 
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Fig. 5.2 Sectional drag coefficient distribution; Advanced Technology Blade. 



Fig. 5.3 Sectional pitching moment coefficient distribution; Advanced Technology Blade. 
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The subscripts (1) and (2) in Eqns. 5.1 - 5.4 correspond to the momentum and force 
equations, respectively, for thrust and torque. This system of equations is then used to 
solve for the inflow and swirl velocities by equating the thrust and torque as follows 


dT i = dT 2 , ( 5 -5) 

dQi= dQ 2 . ( 5 -6) 


The total inflow angle of the blade section (A, in Eqns. 5.2 and 5.4) is defined as 


A — A a + A s } 


(5.7) 


where A a is the angle of the aerodynamic inflow and A s is the additional inflow angle 
which arises due to the inclusion of blade dynamic effects (Fig. 5.4). The effective inflow 
angle due to aerodynamics (A a ) is defined as 


A a 


= tan 


V~ + T>i ) 

Qr-u T j 


(5.8) 


where Voo is the forward velocity of the aircraft and Q. is the rotational speed of the rotor. 



Fig. 5.4 Blade element inflow definitions. 
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To define the effective structural inflow angle, A s , it is first necessary to translate the 
blade dynamics from the global, untwisted coordinate system to a coordinate system which 
is parallel to the effective aerodynamic inflow angle (Fig. 5.5). This is mathematically 
stated as follows. 
v b = v b cos A a + w b sin A a 
w b = - v b sin A a + w b cos A a 

where v b is the dynamic velocity parallel to the aerodynamic inflow angle and w b is the 
velocity perpendicular to the inflow angle. The quantities, v b and w b , represent the 
dynamic effects perpendicular and parallel to the free stream velocity, respectively. These 
transverse velocities are defined as 

v b = vcos(a a ) - wsin(oc a ) 1Q) 

w b = vsin(a a ) + wsin(a a ) - §Y' ac 

where Y ac is the offset between the aerodynamic center of the blade and the axis of twist. 

The physical angle of attack of the blade is denoted a a and the cross-sectional velocities 
in the global, twisted coordinate system are v and w (Fig. 5.5). The effective 
aerodynamic inflow angle is defined as follows. 

_ Voo+^i 

(5-11) 


sin A, 


W 


cosA fl = 


fir - u T 
W 



'N. 



Fig. 5.5 Coordinate systems in blade cross section. 
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where W is the total resultant velocity. Under the assumption that the velocities due to 
airflow are much greater than those due to blade dynamics, the total resultant velocity is 
written as 

W = W a = W s (5-12) 


The effective structural inflow angle is then written as follows. 



(5.13) 


Using Eqns. 5.7 - 5.13, the total inflow angle is expressed as 


A = A a - 


<h( v ~ + v i) , w b (^r-u T ) 

^ r o 


W' 




(5.14) 


The effective angle of attack (a) for the blade cross section is defined as follows 
a = a a + o - A , (5.15) 


where or, is the physical angle of attack of the blade and <j) is the elastic twist due blade 

deformations. Using Eqn. 5.14, the effective angle of attack is rewritten as 

VbO^+Vi) w b (Qr-u T ) 

a= a - A., + 0 + 9 — 7 • 

— 3 W 2 W 2 


Aerodv nanuo 


Elastic deformations 


In the above equation, the first two terms correspond to the purely aerodynamic effects and 
the last three terms represent contributions due to blade elastic deformations. Note that the 
steady state angle (a ss ) is written as 

a ss = a j- -V + 0 . ^ 5 - 17 ) 


5.2 Energy Formulation 

To evaluate the external work due to the aerodynamic loading it is first necessary to 
write the equations for lift, drag and pitching moment, per unit area, as follows. 
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I = |pw 2 % 


|pw 2 c da 


a a “ A a + $ + 




a a -A a +(|) + 


v b (Voo+^i) Wb(^ r ~ u T) > 
W 2 W 2 , 

^(Vc+Di) w b (flr-u T )^ 








J 


m = ^pW 2 cc ma 


a a -A a +(]) + 


v b (Voo +^j) w b (Qr-u T ) 




W‘ 


(5.18a) 

(5.18b) 

(5.18c) 


The coefficients (q a , c da and c ma ) are related to the derivatives of the lift, drag and 

moment coefficients (q, c d , c m ), respectively as follows. 

3c i 


Ci — Ci + 

0 3a 

3c d 

c r ] = c d + — — 

da d ° 3a 


c m a c m 0 + 


3c 


m 


3a 


(5.19a) 

(5.19b) 

(5.19c) 


The external work done by the aerodynamic loads is then written as 
W e = J (-/ w + d v + m (j>) dS 


(5.20) 


where S is the blade surface area and w, v and <() are the global, untwisted displacements. 
From Eqns. 5.18, it must be noted that the aerodynamic forces and moments can be 
separated into a steady aerodynamic term (a a - A a ), a steady term dependent on the elastic 
deflection (<|>) and a term associated with blade dynamics (v b and w b ). Define the 

following parameters 

= --pW 2 (qcosa a + c d sina a ), (5.21a) 


D(j) = — pW 2 (-C] sina a + c d cosa a ), 


(5.21a) 


C C m L(|)Y ac , 


(5.21a) 
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Ha = ~^ pw2 ( Cl a COS0Ca +c dcx sina a)’ 


(5.22a) 


Ha = ^ pW2 Ha Sin0t a + c da cosa a)’ 


M( ^a ~ 9 pW C Cm « Ha Yac ‘ 


(5.22a) 


(5.22a) 


Also, it must be noted that at the outer surfaces of the beam (£i = +h/2) the rotation about 


the x-axis is defined as 


1 ( 3w dv) _ 0w o _ 3w 

2 Oh 3 d ?=h/2 _ an “ 9xi 


(5.23) 


Using Eqns. 5.10 and 5.18, the external work due to steady aerodynamic terms (W ea ) and 


that due to static deflections (W es ) can now be expressed as follows. 


= 


r Ow 

= J(«a- A a) L.jW + D^v + M^— AS, 


(5.24) 


W es=J L 


Ow „ „ 3w Ow , „ 


(5.25) 


The external work due to the dynamics loads (W ed ) is written using Eqns. 5.10, 5.18, 5.20 
and 5.23 as follows. 


W ed -~Jp |c la (vsina + wcosa) + c da (vcosa- wsina) + c c ma 


x Uvcosa- wsina)(Voo + 'Dj)- vsina + wcosa-— (Qr-u T U d 5 


Using the variational principles, the external work functionals are written as follows. 

8W ea = J(a a -A a ) | L <t»-|^( M d} 5 * + D <t>^ d5> 
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where 

a -4 p 


Ci a (V„ + \)i) + c da (Qr-u x ) 


b w = Y' c [c lo ,sina b +c dc( cosa b ](£2r-u x ) 

+ c c ma [(V„ + U;)cosa b - (Or - u x )sina b ] 

c w = pV^[ c i a cosa b -c da sina b ](Qr-u x ) 

- pcc m J(V M +\)j)sina b + (Qr-u x )cosa b ] 


(5.28) 


(5.29) 


(5.30a) 


(5.30b) 


(5.30c) 


5.3 Solution Procedure 

Equations. 5.27 - 5.29 indicate that the external work from the aerodynamic loads will 
yield three separate equations. The first terms are associated with a steady state forcing 
vector independent of the displacements, F aer o (Eqn. 5.27). There are forcing terms that 
are dependent on the displacements (Eqn. 5.28) which will yield a matrix analogous to the 
stiffness matrix ( K aer o)- Finally, there are forcing terms that are dependent on the velocity 
of the displacements (Eqn. 5.29) yielding a damping matrix (C ae ro)- 
Denoting the following matrices 

M *1 T (5.31) 

9q ’ 


and 
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. afv v <tl w w n ] 

11 " dq 


(5.32) 


where q is the degree of freedom vector defined in Chapter 3, the quantities associated with 
the external work due to aerodynamic loads are formulated as follows. 

N 


aero 


i=l 


N 


Jn t f 45 


L3 


Kaero = X 
i=l 


|n t Q] 


d5 


(5.33) 


(5-34) 


and 


N 

c = 'S' 

v "aero / y 
i=l 


Jn 


_ ji 0 0 a w b w 


3-w (^ ^ 


w 


d5 


(5.35) 


where 

f = (a a — A a ) 


3ri 


Da < L, 


K) ! 


Qi = 


3w 


Dv 


- D ^ +2M 


3 2 wl 


*dn 2 


dq 


(5.36) 


(5.37) 


The complete aerodynamic/dynamic/structural equations of motion for the coupled system 
are now written in the following matrix form 

Mq + Cq + Kq = F + (F aero ®^aero *1 ^aero ^)’ (5.38) 


where the quantities within the parentheses correspond to the external work due to 
aerodynamic loading and the remaining terms are associated with the stmctural modeling of 
the beam (Eqn. 3.72) 



106 


5.4 Aeroelastic Stability 

To investigate the aeroelastic stability of the rotor, the natural frequencies of the 
complete coupled equations of motion must be determined. By rearranging the terms in 
Eqn. 5.38, the equation may be rewritten as follows. 

M q + (C - C aero ) q + (K - K aero ) q = F + F aero (5.39) 

The homogeneous portion of Eqn. 5.39 can be rewritten as 

Mq + Cq + Kq = 0, (5-40) 


where 

M = M, 


C 

K 


(C ^aero)’ 

(K-K aero ). 


(5.41) 

(5.42) 

(5.43) 


It is possible to rewrite Eqn. 5.40, which is a second order equation, as a first order 
equation by making the following transformation 



Now Eqn. 5.40 is written as 

{6} = A {5} (5.45) 


where 

r o i 

A = -M -I ic -M -1 C 


(5.46) 


The characteristic roots of Eqn. 5.40 are now determined by solving for the eigenvalues 
of A. It must be noted that the matrix A is neither positive definite, nor symmetric in 
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general and as a result, the eigenvalues are complex. If any of the roots of A have a 
positive real component, the system will be unstable at the given velocity. 

5.5 Correlation 

The adequacy of this aerodynamic representation is demonstrated in Figs. 5.6 and 5.7. 
This relatively simplistic formulation is shown to correlate very well with measured axial 
flow performance of the XV- 15 rotor system in both hover and in airplane mode from tests 
conducted at the Outdoor Aerodynamic Research Facility (O.A.R.F.) and from flight test 
data obtained at NASA Ames Research Center [113,114]. This representation of the rotor, 
which is representative of the original design point of the XV- 15 tiltrotor, is used as the 
baseline, or reference, rotor in this study. Further, as shown in Ref. 13 the results 
obtained using this approach are comparable with those obtained by Dadone et al. [5] in 
which a more comprehensive Euler based analysis technique was used. Therefore, despite 
the relative simplicity, the present approach proves to be quite adequate for modeling prop- 
rotor blade aerodynamics. 






6. Optimization Problem 


The structural, dynamic, aerodynamic and aeroelastic analysis procedures developed 
are now integrated to develop a multidisciplinary optimization procedure for investigating 
the design trade-offs of high speed prop-rotors. The reference aircraft, used as a baseline 
design in the optimization procedure, is a mathematical representation of the XV- 15 tilting 
proprotor aircraft. The rotor is a three-bladed, gimballed rotor with a 25 foot diameter 
[1 13,1 14]. A multipoint optimization procedure is developed and design criteria associated 
with two flight conditions are addresses simultaneously. The first flight condition 
corresponds to hover at sea level and the second flight condition represents high speed 
cruise at an altitude of 25,000 feet and a forward speed of 400 knots. This altitude is 
typical for tiltrotors operating in high speed cruise. 

6.1 Rotor Geometric Modeling 

The rotor planform characteristics are defined as follows. The chord (c) and pre-twist 
angle (0) are defined to have the following cubic spanwise distributions 

C (x) = c 0 + qx + c 2 x 2 + c 3 x 3 , (6.1) 

6(x) = 9j (x - 0. 75) + 0 2 (x - 0. 75) 2 + 0 3 (x - 0. 75) 3 , (6.2) 

where x is the nondimensional radial location (x = x/R, R = blade radius). These cubic 
distributions are selected to give the optimizer sufficient flexibility since the parameters 
which define these distributions are used as design variables. The offset in the twist 
distribution (x - 0.75) is used to ensure zero twist at 75 percent span [117]. The blade 
thickness-to-chord ratio (t/c) is defined to have a quadratic spanwise distribution to ensure a 
monotonic decrease in the thickness from root to tip. This distribution is defined as 
follows 

t/c(x) = to 4- t|X + t2X^. 


(6.3) 
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Similarly, the blade is defined to have a quadratic lifting line offset as follows 

yac<V) = | a oX 2 - ( 6 - 4 ) 

The above distribution is chosen to ensure zero offset at the root as well as zero sweep. 
The blade sweep (A) is defined as follows. 

(6.5) 

or 

A = tan -1 (a 0 x) (6.6) 

6.2 Structural Model 

The box beam dimensions are assumed to be fixed percentages of the chord length and 
airfoil thickness as seen in Fig. 6.1. For the present study, the axis of twist (and rotation) 

is assumed to lie at the centroid of the beam. Further, the beam centroid is assumed to be at 

the 50 percent chord location. 

b(x) = 0.50 c(x) 
d(x) = 0.80 t(x) 



Fig. 6.1 Blade cross section and beam dimensions. 
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The two horizontal walls of the box beam are assumed to have identical composite 
lay-ups. A stacking sequence of {(Pi) 3 /(P2) 3 } 3 is used durin g the optimization. 

Similarly, the two vertical walls are assumed to have the same lay-ups defined as 
1 • A total number of 36 plies is used to ensure that the blade is stiff 


enough to sustain the large aerodynamic loads generated by the rotor. The angle ply 
stacking sequence is selected to investigate the effects of composite ply angles on the 
overall aerodynamic/structural/aeroelastic performance of the rotor. 


6.3 Objective Functions and Constraints 

The optimization problem addresses the simultaneous maximization of both the hover 
figure of merit (FM) and the propulsive efficiency in high speed cruise (r| c ). These 
quantities are defined as follows. 


Pideal 

(6.7) 

~ P 

TVoo 

‘ P 

(6.8) 


The following constraints are imposed to ensure efficient structural and aerodynamic 
performance. To maintain blade aeroelastic stability, constraints are imposed on the real 
part of the stability roots determined from Eqn. 5.39 

X k > -05 k = 1, 2, NAERO (6.9) 


where 05 is the minimum allowable damping, defined to be a small positive number and 
NAERO is the number of modes considered. 

To prevent material failure, constraints are imposed on the individual ply stresses based 
on the Tsai-Wu failure criterion [118]. This criterion assumes that to avoid material failure 
the following equations representing a failure surface in the stress-space must be satisfied. 


Fitfi+FijCjOj < 1 


( 6 . 10 ) 


(i, j = 1, 2, -, 6) 
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where Cj represents the stresses in the coordinate system defined by the material axes (see 
Fig. 6.2). The quantities Fj and Fy are related to tensile and compressive yield strengths of 
the material and are defined as follows. 



( 6 . 11 ) 


( 6 . 12 ) 


where the quantities X, Y are the yield strengths in both compression (subscript *C’) and 
tension (subscript T’) and S is the corresponding shear strength. The quantities F12 and 
F23 are defined as 




(6.14) 


This reduces the total number of constraints as constraints on the individual stress (G\) at 
each ply level are avoided. 
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Fig. 6.2 Composite lamina material axes. 


To maintain rotor thrust at acceptable values during optimization, equality constraints 

are imposed on the thrust in hover and in cruise. These constraints assume the following 

form. 

T h = T href (6-14) 

To = T Cref (6-15) 

6.4 Design Variables 

The design variables that are used during the optimization include the coefficients which 
define the span wise chord and twist distributions (co - C3 and 0i - 03, respectively), the 
thickness-to-chord ratio (t 0 - ti) and the blade sweep (ao). The variables which define the 
composite lay-up in the individual walls (pi - P4) are also used as design variables. It must 
be noted that to ensure realistic blade chord and wing thickness distributions (i.e., positive 
throughout the span), it is necessary to further impose additional geometric constraints on 
these distributions. The minimum allowable nondimensional chord value (c/R) is 
constrained to be 0.02 and the minimum allowable thickness (t) is constrained to be 0.75 
in. Although, the minimum allowable chord values are far too small at the root, due to 
constraints on the stresses these chord distribution is never near critical values except at 
locations near the tip. 



7. Optimization Procedure 


The optimization problem addressed in this research is associated with multidisciplinary 
coupling and involves multiple design objectives and constraints. The Kreisselmeier- 
Steinhauser (K-S) function technique [119] is used to efficiently integrate all of the 
objective functions and constraints into a single envelop function. The problem is thus 
reduced to an unconstrained optimization problem. The Broyden-Fletcher-Goldfarb- 
Shanno (BFGS) algorithm [120] is used to solve the unconstrained nonlinear problem 
(NLP). A hybrid approximate analysis technique based on a two-point exponential 
expansion technique [121] is coupled with the optimization procedure to reduce the 
computational effort. The following sections contain details of the optimization procedure. 

7.1 Kreisselmeier-Steinhauser (K-S) Function Approach 

Since the optimization problem involves more than one design objective, the objective 
function formulation is more complicated. In most of the existing work, the individual 
objective functions are combined using weight factors in a linear fashion [43,44,122]. 
Such methods are judgmental as the answer depends upon the weight factors which are 
often hard to justify. Therefore, the problem is formulated using the Kreisselmeier- 
Steinhauser (K-S) function approach [119]. Using this function the multiple objective 
functions and constraints are transformed into a single envelope function which is then 
minimized using unconstrained optimization techniques. The K-S function has been found 
to perform extremely well by McCarthy et al. in a variety of rotary wing optimization 
problems [9,12-16]. 

The first step in the K-S function approach involves the transformation of the original 
objective functions into reduced objective functions. If the individual objective functions 
are to be minimized, these reduced objective functions assume the following form 

Fk(O) = - 1.0 - gmax ^ 0 . k = 1, - , NOBJ m in (7- la) 
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When the individual objective functions are to be maximized, the reduced objective 
functions are as follows 

Fk(4>) = 1.0 - - gmax < 0 , k = 1, - . NOBJ max (7. lb) 

r Ko 


where Fk 0 represents the value of the original objective function Fk calculated at the 


beginning of each cycle and O is the design variable vector. The quantity gmax is th e value 
of the largest constraint corresponding to the original constraint vector, gj(O) 


(j = 1, 2, ••• , NC) and is held constant during each cycle. These reduced objective 

functions are analogous to constraints, therefore a new constraint vector fmC^) 

(m = 1, 2, , M where M = NC + NOBJ) is introduced which includes the original 

constraints and the constraints introduced through the reduced objective functions (Eqns. 

7.1). The design variable vector in this formulation remains unchanged. The new 

objective function to be minimized is defined using the K-S function as follows 

M 


Fks(O) - tmax + 


jlo 

P " 


(7.2) 


m=l 


where f max is the largest constraint corresponding to the new constraint vector f m (0) and 
in genera! is not equal to g max . The objective function F KS (<E>), which represents an 
envelope function representing the original objective functions and constraints, can now be 
minimized using an\ unconstrained optimization technique. 

The optimization algorithm, based upon this technique, can be explained as follows. 
Initially in an infeasible design space, where the original constraints are violated, the 
constraints due to the reduced objective functions (Eqns. 7.1) are satisfied, i.e. gmax is 
negative. Once the original constraints are satisfied, the constraints due to the reduced 
objective functions become violated. When this happens, the optimizer attempts to satisfy 
these constraints and in doing so, minimizes (or maximizes) the original objective functions 
(Fk). The multiplier p in Eqn. 7.2 is analogous to a draw-down factor where p controls 
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the distance from the surface of the K-S objective function to the surface of the maximum 
constraint function. When p is large, the K-S function closely follows the surface of the 
largest constraint function and when p is small, the K-S function includes contributions 
from all violated constraints. Additional details can be found in Refs. 93, 94, 97-101. 

The K-S function formulation is illustrated for a problem where a single objective 
function is to be maximized subject to two constraints using one design variable (Figs. 7.1 
and 7.2). An initial design point of X = 7 is used in the example. At this point, both 
constraints are satisfied and gmax is therefore negative. As a result, the reduced objective 
function is positive and the constraint introduced through Eqn. 7.1b is violated (Fig. 7.1). 
The three constraints of the problem introduced by the original constraints and the reduced 
objective function are shown in Fig. 7.2 along with the associated K-S function for three 
different values of p. As seen from the figure, for a value of p =1, the K-S function 
represents a more composite envelope function which includes contributions from all three 
constraints. This is especially evident at locations where the values of two or more 
constraints are very similar. For the larger values of p = 3 and p = 5, the K-S function 
envelope more closely represents only the largest constraint even at locations where the 
constraints are similar in value. This simple example demonstrates how larger values of p 
“draw down” the K-S function closer to the value of the largest constraint. 

The K-S function (Fig. 7.2) is minimized using standard nonlinear, unconstrained 
optimization techniques. Once a local minima is reached, a new cycle begins with the 
calculation of a new value of g max and the reformulation of the reduced objective functions 
and the K-S function. The process is repeated until either the original objective functions 
or the design variable vector converges. Note that it is permissible to allow the value of p 
to change from cycle to cycle. This is typically done is a monotonically increasing manner 
so that as the optimization proceeds, the K-S function more closely represents only the 
largest constraint (or the reduced objective function). 



constraint 1 



Fig. 7.1 K-S function formulation example, 
constraint 1 constraint 2 
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7.2 Approximate Analysis 

The two-point exponential approximation technique [121], which was found to perform 
well in other nonlinear optimization problems [8-12, 14-16], is used for the approximation 
of the objective functions and the constraints. This technique derives its name from the fact 
that the exponent used in the expansion is based upon gradient information from the 
previous and current design cycles. The technique is described below. 


£ k (<t>) = F k (<D 0 ) + 



- 1.0 


*°. 3F k (O 0 ) 

Pn 


(7.3) 


where £k(<b) is the approximation of the original objective function Fk(0 0 ). The 
approximate values for the constraints, gj(3>), are similarly calculated. The exponent, p n , 
is defined below 


>°gej 

[aF(4>,)l 

l J 

>->OgJ 

f9F(<S> 0 )l 

l J 

l°g e | 

fo. 1 

L A n J 


c£> ! 

L °n J 



(7.4) 


where the quantity <&i refers to the design variable vector from the previous iteration and 
the quantity <f> 0 denotes the current design vector. A similar expression is derived for the 
constraint vector. The exponent p n can be considered as a “goodness of fit” parameter, 
which explicitly determines the trade-offs between traditional and reciprocal Taylor series 
based expansions. Therefore, the procedure can also be regarded as a hybrid 
approximation technique. It can be seen from Eqn. 7.3 that in the limiting case of p n = 1, 
the expansion is identical to the traditional first order Taylor series and when p n = -1, the 
two-point exponential approximation reduces to the reciprocal expansion form. The 
exponent is then defined to lie within this interval. Therefore, if the exponent p n > 1, it is 
set identically equal to one and if p n < -1, it is set equal to -1. From Eqns. 7.3 and 7.4, it 
is obvious that many singularity points exist in the use of this method. Therefore, care 
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must be taken to avoid such points. In the present study, the linear Taylor series 
approximation is used at such singular points. 

To ensure the validity of the approximation it is necessary to impose bounds, or “move 
limits” on the design variables during the optimization so that the design point remains in 
the neighborhood of the original point. These move limits represent a percent change from 
the original design variable. The move limits in this study are calculated based on a 
variable scheme developed by Thomas et al. [123]. This algorithm adjusts the values of the 
move limits based on changes in the maximum violated constraint and also by tracking the 
individual move limits to see whether they reach the same upper or lower limit over three 
consecutive evaluations. Another important aspect of the scheme developed in Ref. 123 is 
the ability to allow design variables to cross over between negative and positive values. 
Details of move limit approach are presented in Ref. 123. 



8. Optimization Results 


The reference rotor used is representation of the XV- 15 proprotor which is an advanced 
three-bladed gimballed rotor [1 13,1 14]. The aerodynamic optimization is performed at a 
cruise altitude of 25,000 feet and a forward velocity of 400 knots with a rotational speed of 
421 RPM. A vehicle weight of 13,000 lbs and aircraft lift to drag ratio (L/D) of 8.4 is 
assumed. Therefore, the thrust in cruise is constrained to be at 774 lbs for the two engine 
aircraft. In hover, the aircraft is assumed to be operating at sea level conditions with a 
rotational speed of 570 RPM and a 12 percent down load effect from the rotor/wing 
interaction. The thrust in hover is therefore constrained to be at 7280 lbs. The blade is 
discretized into 10 aerodynamic segments (11 node points) and the composite box beam is 
similarly discretized into 10 spanwise elements and 1 chordwise element for a total of 564 
degrees of freedom. The composite material used in the structural analysis is carbon-PEEK 
AS4/APC2 [124] and the material properties are presented in Table 8.1. A total of 15 
design variables are used during the optimization. The optimum design converges in 35 
cycles and the results are presented in Table 8.2 and Figs. 8.1 - 8.6. 

Table 8.1 Summary of beam material properties 

Ei = 19.4 x 10 6 p.s.i., E 2 = 1.29 x 10 6 p.s.i., 

G 12 = 0.740 x 10 6 p.s.i., Gb = 0.500 x 10 6 p.s.i., 
p 12 = 0.28, p = 1.80 x 10-3 slug/in 3 

ply thickness = 0.001 in. 

Ultimate Strengths 

X T = 309 x 10 3 p.s.i., X c = 160 x 10 3 p.s.i., 

Y t = 11.6 x 10 3 p.s.i., Y c = 29.0 x 10 3 p.s.i., 

S = 23.2 x 10 3 p.s.i. 
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From Table 8.2 and Fig. 8.1 it is seen that the hover figure of merit (FM) increases by 
3.6 percent and the high speed cruise propulsive efficiency (r|ax) is significantly improved 
(55 percent). It must be noted, however, that the baseline rotor was originally designed for 
operation at 300 knots and therefore has a poor cruise efficiency (p a x = 0.49) at the 
optimization design speed of 400 knots. As a result, the improvement in Pax is much more 
significant compared to the increase in the hover figure of merit which has a fairly high 
value initially. A complete understanding of the aerodynamic improvements is obtained by 
examining the design variable trends. 

Table 8.2 Comparison of optimum results 

Reference Optimum 

Objective Functions 


FM 

0.7691 

0.7974 

tlax 

0.4856 

0.7502 

Design Variables 

c 0 

0.1094 

0.1050 

ci 

-0.09256 

-0.09760 

C2 

0.1575 

0.1630 

c 3 

-0.08176 

-0.07630 

01 (rad.) 

-0.3455 

-0.3443 

02 (rad.) 

0.7693 

0.5817 

03 (rad.) 

-0.1461 

-0.1057 

to 

0.3155 

0.2167 

ti 

-0.3193 

-0.2370 

t2 

0.07517 

0.09128 

ao 

0.000 

-0.05615 

Pi (deg.) 

0.0 

0.0 

p2 (deg.) 

0.0 

0.5 

p3 (deg.) 

0.0 

-1.7 

P 4 (deg.) 

0.0 

-1.5 
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□ Reference 



Optimum 



Fig. 8.1 Comparison of optimum results. 


1 


The optimum chord distribution is presented in Fig. 8.2 and shows that the chord 
distribution over the majority of the blade span is reduced from the baseline values. At the 
tip, however, the chord is actually increased. It is important to note that although the area 
weighted solidity io A ) is slightly decreased from the reference value (2.8 percent), both the 
thrust weighted solidity (g t ) and the power weight solidity (a P ) are not decreased as much 
(2.3 percent and 1.8 percent, respectively) as shown in Fig. 8.3. The solidity ratios are 
defined as follows. 

<5:=^-, i = A, T or P, (8.1) 

1 u R 


where 
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( 8 . 2 ) 


(8.3) 


(8.4) 


In the above equations, c is blade chord and x is the nondimensional radius. The reason 
for the reduced chord near the root is due to the fact that in the reference blade this section 
produces a significant amount of drag in cruise (Fig. 8.4) without generating any 
significant lift (Fig. 8.5). After optimization, it is seen that the drag is significantly reduced 
in this region whereas the lift is only slightly affected. In hover, the root section generates 
very little lift and drag (Figs. 8.6 and 8.7, respectively) therefore the root chord is reduced 
only due to constraints on the stresses. In the absence of the stress constraints, it is likely 
that the root chord would be reduced to smaller values. As a result, the optimizer reduces 
the chord throughout the blade span, except near the tip. It is interesting to note, that in the 
optimum configuration, there is an increase in the chord values from about midspan 
towards the tip, resulting in a slight inverse taper. This is due to the fact that this outboard 
section represents the working section of the blade. In both hover and in cruise, the 
majority of the lift is produced by the outer 25 percent of the blade. By increasing the 



Rotor solidity 
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chord in this region of the blade, the optimizer is redistributing the load to a region which is 
beneficial to both flight conditions. 



Fig. 8.2 Comparison of chord distributions. 
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Fig. 8.3 Comparison of rotor solidity. 
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Fig. 8.4 Comparison of high speed cruise sectional drag. 



Fig. 8.5 Comparison of high speed cruise sectional lift. 
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Nondimensional radius, x/R 


Fig. 8.6 Comparison of hover sectional lift. 
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Fig. 8.7 Comparison of hover sectional drag. 
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The optimum and the reference blade thickness distributions are presented in Fig. 8.8. 
From the figure it is seen the blade thickness is reduced from the reference values 
throughout the entire span, except at the tip. One reason for this reduction is to reduce the 
profile drag of the blade by reducing the thickness. Reductions in the profile drag, in turn 
improves the aerodynamic performance. A second reason for this reduction is to increase 
the drag divergence Mach number (Mdd). This Mach number is defined as the point where 
a further increase in Mach number will result in a sharp increase in the drag [111]. In the 
reference blade the local Mach numbers near the tip in cruise are very near Mdd- Through 
reductions in the blade thickness distributions, the optimizer increases Mdd throughout the 
blade and as a result the drag is reduced improving the performance of the rotor. The 
constraints imposed on the blade thickness at the tip prevent the optimizer from reducing 
the tip thickness below the reference value. Near the tip, these constraints become active 
and as a result no further reductions in blade thickness are obtained in this region. 



Fig. 8.8 Comparison of thickness-to-chord distributions. 
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Figure 8.9 presents the twist distributions of both the reference and the optimum rotors 
where it is seen that the twist is reduced from root to about midspan of the rotor. These 
reductions are due to the fact that the optimizer is attempting to unload this section of the 
blade since the drag in the cruise condition is very high for the reference blade (Fig. 8.4). 
The effective angle of attack is lower as a result of the reduction in blade twist and this 
results in reductions in both lift and drag in the region. This is again due to the fact that the 
optimizer is redistributing the lift outboard towards the tip to improve the performance. 
Near the tip, the twist distributions are very similar. 



Fig. 8.9 Comparison of twist distributions. 

The lifting line offset (Y ac ) and the resulting sweep distribution ( A ) are presented in 
Figs. 8.10 and 8.11. These figures indicate that very little sweep is introduced after 
optimization and that the blade is swept forward. This can be explained as follows. The 
introduction of sweep reduces the effective Mach number which in turn improves the 
performance of the high speed cruise propulsive efficiency. However, only a slight 
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amount of sweep is introduced (about 3 degrees at the tip) and the reductions in the 
effective Mach numbers are not significant. This is due to the fact that in cruise there is a 
large nose down moment throughout the blade span. These moments are larger than the 
corresponding moments introduced through the lifting distribution and the offset of the 
aerodynamic center from the axis of twist. As a result, the blade twists down in cruise. 
However, in hover, the moments due to the lifting distribution and the aerodynamic offset 
are much larger than the corresponding moments due to the aerodynamic pitching 
moments. Therefore the blade twists backward (nose up) in hover. By sweeping the blade 
slightly forward, the offset between the lift and the axis of twist is increased and this 
reduces the amount of negative twist in cruise which in turn improves the performance. 
However, the forward sweep increases the amount of positive twist in hover and as a result 
only a slight amount of twist can be allowed without adversely affecting the hover 
performance. Further, a significant amount of sweep would increase the bending moments 
which would result in increased blade stresses. 



Fig. 8.10 Comparison of lifting line distributions. 
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Fig. 8.11 Comparison of blade sweep distributions. 

The stability roots of the five most critical modes are presented in Figs. 8.12 (hover) 
and 8.13 (cruise). The figures show that in both flight conditions the real part of the 
stability root is negative for both the reference and optimum rotors assuring that the system 
is stable. It is further seen in the figures that an imaginary line drawn through the locus of 
roots, originating at the origin, is nearly linear. This is due to the fact that the damping due 
to aerodynamics is small. Therefore, the structural damping (assumed to be two percent, 
proportional damping) dominates the aerodynamic damping. This is caused by several 
factors. First, due to the large loads under which the blade operates (each blade must 
generate nearly 2500 lbs of thrust), the blade must be extremely stiff in order to withstand 
the bending stresses. Also, unlike fixed wing aircraft, the aerodynamic loading is not 
significantly altered with changes in forward speed due to the fact that the rotor is trimmed 


to a constant thrust value. 
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Fig. 8.12 Comparison of hover stability roots; 2 percent structural damping. 
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Fig. 8.13 Comparison of high speed cruise stability roots; without structural damping. 
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The optimum ply angles are presented in Table 8.2. Very little changes are observed 
between the reference and the optimum values. This is due to the fact that the stability 
constraints are never critical during optimization. Since the ply angles have very little effect 
on the aerodynamic performance and the stability roots are never near critical, there is little 
change in the ply angles after optimization. Figures 8.14 - 8.21 present the midplane 
stresses in the individual walls, at each comer, before and after optimization. These figures 
illustrate the increase in the root stresses in each wall for both flight conditions and show 
that these stresses are more critical in hover (Figs. 8.14 - 8.17) than in cruise (Figs. 8.18 - 
8.21) . The increase in the stresses once again due to the reduction in cross-sectional area 
which is caused by the reductions in root chord. Although, the stresses are increased after 
optimization, the overall Tsai-Wu criterion is satisfied at each ply level at every comer for 
both flight conditions. In the absence of these stress constraints, it is likely that the 
optimizer would have further reduced the root chord in order to improve the aerodynamic 
performance. 

To fully investigate the phenomenon of tilt-rotor aeroelastic stability, the combined 
problem of the wing/rotor/pylon assembly of the tilt-rotor aircraft should be investigated. 
Further, the use of a more comprehensive aerodynamic analysis, such as panel code, may 
be beneficial. If such measures were taken, it is believed the ply angles (particularly in the 
aircraft wing) would have played a more significant role on the overall performance and 
stability of tilt-rotors. 
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Fig. 8.16 Comparison of midplane stresses; bottom wall, hover. 
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Fig. 8.17 Comparison of midplane stresses; left wall, hover. 
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Fig. 8.18 Comparison of midplane stresses; top wall, high speed cruise. 
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Fig. 8.19 Comparison of midplane stresses; right wall, high speed cruise. 





Fig. 8.20 Comparison of midplane stresses; bottom wall, high speed cruise. 



Fig. 8.21 Comparison of midplane stresses; left wall, high speed cruise. 



9. Concluding Remarks 


A new beam theory has been developed to model composite box beams with arbitrary 
wall thicknesses. The theory, which is based on higher-order composite laminate theory, 
approximates the three-dimensional elasticity solution rather than reducing the cross- 
sectional properties to one-dimensional beam properties. Arbitrary spanwise distributions 
of blade twist, taper and sweep are included in the formulation. The developed theory 
satisfies the stress free boundary conditions on the inner and outer surfaces of the beam. 
Both inplane and out-of-plane warping are included in the formulation. 

Next a procedure for the aeroelastic stability of prop-rotor blades has been developed. 
The aerodynamic loads are based on the classical blade element momentum theory. The 
coupled equations of motion are developed which represent a trimmed blade configuration. 
Finally, the developed structural, aerodynamic and aeroelastic procedure are integrated 
within an optimization procedure to investigate prop-rotor performance in both high speed 
cruise and hover. The optimization problem includes multiple objectives and the 
Kreisselmeier-Steinhauser (K-S) function is used to formulate the optimization problem. 
This function represents an envelope function of all of the objective functions and 
constraints and transforms the problem to an unconstrained optimization problem. The 
Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is used to determine the search 
direction. The procedure is coupled with a hybrid approximate expansion is used to reduce 
the computational effort. The following important observations are made. 

1 . Very good overall agreement is observed between the static results and available 
experimental data for thin-walled beams. The dynamic results correlate well with 
NASTRAN using both solid elements and shell elements. 

2 . For large angle ply laminates (e.g. 45°), the present approach predicts the behavior 
very well as shown by the correlation with experimental results. For these 
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laminates, the present approach represents a significant improvement from a 
previously developed quasi-analytical method. 

3 . The results from the symmetric beams have better overall correlation than the anti- 
symmetric beams. This is due to the fact that magnitude of the twist measured for 
the anti-symmetric beam correlation is very small and therefore is more difficult to 
determine experimentally. 

4. The effect of transverse shear stresses is critical in case of thick-walled sections. 
This introduces large nonlinearities in the twist distribution. Further, the local twist 
in the individual walls is not equal as it is in case of thin- walled beams. 

5 . The “beam” twist is a local quantity which can only be defined at a point in the 
cross section. Arbitrary definition of the twist at a convenient point in the beam 
cross-section is inaccurate for thick-walled cross sections. 

6. The modes shapes often display a significant amount of bending-twist coupling 
and/or extension-shear coupling. The coupling is more noticeable in beams with 
thicker wall sections. 

7 . The present beam theory captures the effects of inplane and out-of-plane warping. 
For thick-walled beams with low aspect ratios, the warping terms are significant 
even at the lower modes. For the thin-walled beams, the inplane warping is more 
important than the out-of-plane warping. Results obtained using NASTRAN tend 
not to capture the lower warping modes. 

8. The increased warping in beams with thicker walls is due to the presence of 
transverse shear stresses through the thickness of the walls which increases with 


laminate thickness. 
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9 . The multidisciplinary optimization developed yields significant improvements in the 
aerodynamic performance without any degradation of the structural response and 
aeroelastic stability. 

10. Improvements in the aerodynamic performance are obtained largely through 
changes in the blade chord, twist and thickness distribution. By changing the chord 
and twist, the load is redistributed to the outboard sections of the blade. Reductions 
in the blade thickness reduce the profile drag thereby increasing the performance. 

1 1 . The optimizer does not significantly alter the blade sweep due to the conflicting 
requirements between the two flight conditions. 

12. Changes in the ply angles are negligible due to the fact that aeroelastic stability 
constraints are never critical to optimization. 

13. Using the current aerodynamic analysis, the aeroelastic stability of prop-rotor 
blades is not a significant issue due to the fact blade must be extremely stiff in order 
to withstand bending stresses and the fact that the blade is trimmed to a constant 
thrust value. A more complete representation of the aircraft, including the 
wing/rotor/pylon assembly is recommended in order to better evaluate the 
aeroelastic stability of tilt-rotor aircraft. 
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The transformation between the untwisted global coordinate system and the rotated 


coordinate system is expressed as follows. 
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The two unswept wall coordinate systems are written in terms of the two beam coordinate 


systems as 
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where the global coordinate system (Xi, Yi, Z\) is written with the subscript T to indicate 
the fact that Yj is always aligned with the reference (untwisted) width of the individual 
walls. Similarly, the rotated coordinate system (X- , Y-, Z-) is always aligned such that Y- 


is parallel to the rotated width of the individual wall. This notation is adopted so that the 

relationship between the global and local coordinate systems can be expressed using only a 
single set of equations. The quantities, Y oi , Z oi , Y'. and Z^. correspond to distances 


from the axis of rotation to the edge of the individual walls within which the wall 
coordinate systems are defined. Finally, the relationship between the local, swept and 
twisted coordinate system and the local, unswept and twisted coordinate system is written 
as 
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where z 0 . is the offset between the swept and unswept coordinate systems. Note that in 

the horizontal walls (i = 1 and 3), this offset is zero. Using Eqns. A1 - A4, the local 
wall, twisted coordinate system is written as 
1 

ii 1 i i i 
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It must be noted that in the above relationship, both 6 and z 0 . are functions of x and vary 

along the span. The Jacobian matrix between the local untwisted and local twisted 
coordinate systems is then expressed as 
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where the subscript T has been omitted for convenience. It is seen from Eqn. A6, that the 
determinant of the Jacobian is equal to one. 
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The zeroth order inplane strains in the absence of pre-twist are defined as follows. 
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and the third components of the inplane strain are 
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and the second order components of the out-of-plane strains are defined as 
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The additional non-zero inplane strain components due to pre-twist are written as 
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The non-zero out-of-plane components 
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of the strain due to presence of pre-twist are written 
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Finally, the nonzero strain terms associated with sweep (which exist only in the vertical 


walls) are written as follows. 
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Next, the convergence study is conducted in the natural frequencies and a 
nondimensional natural frequency parameter, is used where 


= COj 



(C.l) 


D = Eh 3 /l2(l- v 2 ). 


(C.2) 


The natural frequency parameters, for the first six modes are presented for a range of 
mesh sizes in Figs C.2 - C.7. From these figures it is again noted that although 
convergence can be achieved with only one or two chordwise elements, often the 
converged results are significantly different from the actual results. This is particularly 
evident in the third flapping mode (Fig. C.5), the second torsional mode (Fig. C.6) and the 
first plate mode (Fig. C.l). For a mesh with three or more chordwise elements, 
convergence is achieved very quickly. In general for such chordwise discretizations, 
convergence is achieved with 10-12 spanwise elements. 



Fig. C.l Bending slope for 30° twisted, isotropic plate subjected to 1 lb. tip bending load. 
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Natural frequency parameter for the first flapping mode; 30° twisted, isotropic 

plate. 



Fig. C.3 Natural frequency parameter for the second flapping mode; 30° twisted, isotropic 

plate. 
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Fig. C.4 Natural frequency parameter for the first torsion mode; 30° twisted, isotropic 

plate. 



Fig. C.5 Natural frequency parameter for the third flapping mode; 30° twisted, isotropic 

plate. 
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Fig. C.6 Natural frequency parameter for the second torsion mode; 30° twisted, isotropic 

plate. 



Fig. C.7 Natural frequency parameter for the first plate mode; 30° twisted, isotropic plate. 
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C.2 Composite, 45° symmetric beam 

Convergence tests for the beam are performed for both a static case with a 1 lb. tip 
bending load as well as for the first several natural frequencies. The beam used for this 
study is a symmetric 45° graphite/epoxy composite beam for which the properties are listed 
in Table. C.2. Figures C.8 - C. 10 present the static deflections at the tip. It must be noted 
that the values presented in these figures represent the average value for the entire cross 
section. Since the developed theory does not rely on reducing the beam behavior to one- 
dimensional parameters, values of the cross-sectional flap, lag and twist are not explicitly 
calculated. However, in case of thin-walled beams, the values in each of the walls are very 
close to each other, therefore, in such cases the approximation of the cross-sectional twist 
being equal to average twist of all four walls is valid. It is observed in Figs. C.8 - C.10 
that convergence is achieved with approximately 10 spanwise elements. More importantly 
it is noted that this convergence is independent of the number of chordwise elements. This 
is largely due to the fact that the beam is decomposed into the individual walls which make 
up the beam and as a result, a single chordwise element in each wall results in four 
chordwise element for the beam cross section. 

Table C.2 Details of symmetric 45° composite beams 

Flanges ” Webs 

Top Bottom Left Right 

[4 ~6 W\~6 [457-15°] 3 [457-45 °] 3 

Material Properties 

E l = 20.59 x 10 6 p.s.i., E t = 1.42 x 10 6 p.s.i., 

G L t = 0.89 x 10 6 p.s.i., Vlt = 0.42. 

Beam Dimensions 

Length = 30 in., width = 0.953 in., depth = 0.53 in., 
ply thickness = 0.005 in, number of plies = 6, 

wall thickness = 0.030 in. 
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Similar trends are observed in the convergence rates of the natural frequencies which 
are presented in Figs. C.ll - C.17. It is seen from these figures that once again 
convergence is independent of the number of chordwise elements. For the lower 
frequencies, convergence is obtained with approximately 10 spanwise elements. For the 
higher frequencies, convergence is a little slower and although the results are not truly 
converged with 10 spanwise elements, they are quite close to the converged values. 



Fig. C.8 Flapping displacement for 45° symmetric beam subjected to 1 lb. tip bending 

load. 
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Fig. C.9 Lagging displacement for 45° symmetric beam subjected to 1 lb. tip bending load. 



Fig. C.10 Elastic twist for 45° symmetric beam subjected to 1 lb. tip bending load. 
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